Changeset d368d21 in sasmodels
- Timestamp:
- Feb 6, 2016 10:56:54 AM (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:
- 321736f
- Parents:
- 03582f9 (diff), d5e650d (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Files:
-
- 34 added
- 52 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/guide/index.rst
r19dcb933 rbb6f0f3 1 ********************** 2 SAS Model Organization 3 ********************** 1 ********** 2 SAS Models 3 ********** 4 5 Small angle X-ray and Neutron (SAXS and SANS) scattering examines the 6 scattering patterns produced by a beam travelling through the sample 7 and scattering at low angles. The scattering is computed as a function 8 of $q_x$ and $q_y$, which for a given beam wavelength corresponds to 9 particular scattering angles. Each pixel on the detector corresponds to 10 a different scattering angle. If the sample is unoriented, the scattering 11 pattern will appear as rings on the detector. In this case, a circular 12 average can be taken with 1-dimension data at $q = \surd (q_x^2 + q_y^2)$ 13 compared to the orientationally averaged SAS scattering pattern. 4 14 5 15 Models have certain features in common. -
extra/pylint.rc
r5ef0633 r823e620 6 6 # Python code to execute, usually for sys.path manipulation such as 7 7 # pygtk.require(). 8 #init-hook= 8 #init-hook='import os, sys; sys.path.extend([os.getcwd(),os.path.join(os.getcwd(),"extra")])' 9 init-hook='import sys; sys.path.append('extra') 9 10 10 11 # Profiled execution. 11 profile=no12 #profile=no 12 13 13 14 # Add files or directories to the blacklist. They should be base names, not … … 64 65 star-args, 65 66 unbalanced-tuple-unpacking, 67 locally-enabled, 66 68 locally-disabled, 67 69 … … 90 92 # Add a comment according to your evaluation note. This is used by the global 91 93 # evaluation report (RP0004). 92 comment=no94 #comment=no 93 95 94 96 # Template used to display messages. This is a python new-style format string … … 99 101 100 102 # Required attributes for module, separated by a comma 101 required-attributes=103 #required-attributes= 102 104 103 105 # List of builtins function names that should not be used, separated by a comma … … 106 108 # Good variable names which should always be accepted, separated by a comma 107 109 good-names=_, 108 input, 109 i,j,k,n,x,y,z, 110 q,qx,qy,qz, 111 dt,dx,dy,dz,id, 112 Iq,dIq,Qx,Qy,Qz, 113 p, 110 input, id, 111 h, i, j, k, n, p, v, w, x, y, z, 112 q, qx, qy, qz, Q, Qx, Qy, Qz, 113 dt, dx, dy, dz, dq, 114 Iq, dIq, Iqxy, Iq_calc, 114 115 ER, call_ER, VR, call_VR, 116 Rmax, SElength, 115 117 116 118 # Bad variable names which should always be refused, separated by a comma … … 278 280 # (useful for modules/projects where namespaces are manipulated during runtime 279 281 # and thus existing member attributes cannot be deduced by static analysis 280 ignored-modules=numpy,np 282 ignored-modules=numpy,np,numpy.random, 283 bumps,sas, 281 284 282 285 # List of classes names for which member attributes should not be checked … … 286 289 # When zope mode is activated, add a predefined set of Zope acquired attributes 287 290 # to generated-members. 288 zope=no291 #zope=no 289 292 290 293 # List of members which are set dynamically and missed by pylint inference … … 316 319 # List of interface methods to ignore, separated by a comma. This is used for 317 320 # instance to not check methods defines in Zope's Interface base class. 318 ignore-iface-methods=isImplementedBy,deferred,extends,names,namesAndDescriptions,queryDescriptionFor,getBases,getDescriptionFor,getDoc,getName,getTaggedValue,getTaggedValueTags,isEqualOrExtendedBy,setTaggedValue,isImplementedByInstancesOf,adaptWith,is_implemented_by321 #ignore-iface-methods=isImplementedBy,deferred,extends,names,namesAndDescriptions,queryDescriptionFor,getBases,getDescriptionFor,getDoc,getName,getTaggedValue,getTaggedValueTags,isEqualOrExtendedBy,setTaggedValue,isImplementedByInstancesOf,adaptWith,is_implemented_by 319 322 320 323 # List of method names used to declare (i.e. assign) instance attributes. … … 335 338 336 339 # Maximum number of arguments for function / method 337 max-args= 5340 max-args=15 338 341 339 342 # Argument names that match this expression will be ignored. Default to name … … 342 345 343 346 # Maximum number of locals for function / method body 344 max-locals= 15347 max-locals=25 345 348 346 349 # Maximum number of return / yield for function / method body -
extra/pylint_numpy.py
r63b32bb r823e620 8 8 #print("processing",module.name) 9 9 if module.name.startswith('numpy'): 10 if module.name == 'numpy': import numpy 11 elif module.name == 'numpy.random': import numpy.random 10 if module.name == 'numpy': 11 import numpy 12 elif module.name == 'numpy.random': 13 import numpy.random 14 from numpy.random import seed, get_state, set_state 12 15 -
extra/pylint_pyopencl.py
r3c56da87 r823e620 9 9 if module.name == 'pyopencl': 10 10 import pyopencl 11 import pyopencl as cl -
extra/run-pylint.py
rf01e53d r823e620 6 6 7 7 def main(): 8 extra = abspath(dirname(__file__)) 9 root = abspath(joinpath(extra, '..')) 10 8 11 envpath = os.environ.get('PYTHONPATH',None) 9 12 path = [envpath] if envpath else [] 10 path.append(abspath(dirname(__file__))) # so we can find the plugins 13 path.append(extra) 14 15 #bumps = abspath(joinpath(root, "..", "bumps")) 16 #periodictable = abspath(joinpath(root, "..", "periodictable")) 17 #sasview = abspath(joinpath(root, "..", "sasview", "src")) 18 #path.extend((bumps, periodictable, sasview)) 19 11 20 os.environ['PYTHONPATH'] = ':'.join(path) 12 root = abspath(joinpath(dirname(__file__), '..')) 21 22 # Run the lint command 23 cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels" 13 24 os.chdir(root) 14 cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels"15 25 status = os.system(cmd) 16 26 sys.exit(status) -
sasmodels/__init__.py
r9890053 r37a7252 1 """ 2 sasmodels 3 ========= 4 5 **sasmodels** is a package containing models for small angle neutron and xray 6 scattering. Models supported are the one dimensional circular average and 7 two dimensional oriented patterns. As well as the form factor calculations 8 for the individual shapes **sasmodels** also provides automatic shape 9 polydispersity, angular dispersion and resolution convolution. SESANS 10 patterns can be computed for any model. 11 12 Models can be written in python or in C. C models can run on the GPU if 13 OpenCL drivers are available. See :mod:`generate` for details on 14 defining new models. 15 """ 16 1 17 __version__ = "0.9" -
sasmodels/bumps_model.py
rec7e360 r190fc2b 6 6 arguments to set the initial value for each parameter. 7 7 8 :class:`Experiment` combines the *Model* function with a data file loaded by the9 sasview data loader. *Experiment* takes a *cutoff* parameter controlling8 :class:`Experiment` combines the *Model* function with a data file loaded by 9 the sasview data loader. *Experiment* takes a *cutoff* parameter controlling 10 10 how far the polydispersity integral extends. 11 11 … … 19 19 from .direct_model import DataMixin 20 20 21 __all__ = [ 22 "Model", "Experiment", 23 ] 24 21 25 # CRUFT: old style bumps wrapper which doesn't separate data and model 26 # pylint: disable=invalid-name 22 27 def BumpsModel(data, model, cutoff=1e-5, **kw): 28 r""" 29 Bind a model to data, along with a polydispersity cutoff. 30 31 *data* is a :class:`data.Data1D`, :class:`data.Data2D` or 32 :class:`data.Sesans` object. Use :func:`data.empty_data1D` or 33 :func:`data.empty_data2D` to define $q, \Delta q$ calculation 34 points for displaying the SANS curve when there is no measured data. 35 36 *model* is a runnable module as returned from :func:`core.load_model`. 37 38 *cutoff* is the polydispersity weight cutoff. 39 40 Any additional *key=value* pairs are model dependent parameters. 41 42 Returns an :class:`Experiment` object. 43 44 Note that the usual Bumps semantics is not fully supported, since 45 assigning *M.name = parameter* on the returned experiment object 46 does not set that parameter in the model. Range setting will still 47 work as expected though. 48 49 .. deprecated:: 0.1 50 Use :class:`Experiment` instead. 51 """ 23 52 warnings.warn("Use of BumpsModel is deprecated. Use bumps_model.Experiment instead.") 53 54 # Create the model and experiment 24 55 model = Model(model, **kw) 25 56 experiment = Experiment(data=data, model=model, cutoff=cutoff) 26 for k in model._parameter_names: 27 setattr(experiment, k, getattr(model, k)) 57 58 # Copy the model parameters up to the experiment object. 59 for k, v in model.parameters().items(): 60 setattr(experiment, k, v) 28 61 return experiment 29 62 63 30 64 def create_parameters(model_info, **kwargs): 65 """ 66 Generate Bumps parameters from the model info. 67 68 *model_info* is returned from :func:`generate.model_info` on the 69 model definition module. 70 71 Any additional *key=value* pairs are initial values for the parameters 72 to the models. Uninitialized parameters will use the model default 73 value. 74 75 Returns a dictionary of *{name: Parameter}* containing the bumps 76 parameters for each model parameter, and a dictionary of 77 *{name: str}* containing the polydispersity distribution types. 78 """ 31 79 # lazy import; this allows the doc builder and nosetests to run even 32 80 # when bumps is not on the path. 33 81 from bumps.names import Parameter 34 35 partype = model_info['partype']36 82 37 83 pars = {} … … 40 86 value = kwargs.pop(name, default) 41 87 pars[name] = Parameter.default(value, name=name, limits=limits) 42 for name in partype['pd-2d']:88 for name in model_info['partype']['pd-2d']: 43 89 for xpart, xdefault, xlimits in [ 44 ('_pd', 0., limits),45 ('_pd_n', 35., (0, 1000)),46 ('_pd_nsigma', 3., (0, 10)),47 ]:90 ('_pd', 0., limits), 91 ('_pd_n', 35., (0, 1000)), 92 ('_pd_nsigma', 3., (0, 10)), 93 ]: 48 94 xname = name + xpart 49 95 xvalue = kwargs.pop(xname, xdefault) … … 51 97 52 98 pd_types = {} 53 for name in partype['pd-2d']:99 for name in model_info['partype']['pd-2d']: 54 100 xname = name + '_pd_type' 55 101 xvalue = kwargs.pop(xname, 'gaussian') … … 63 109 64 110 class Model(object): 111 """ 112 Bumps wrapper for a SAS model. 113 114 *model* is a runnable module as returned from :func:`core.load_model`. 115 116 *cutoff* is the polydispersity weight cutoff. 117 118 Any additional *key=value* pairs are model dependent parameters. 119 """ 65 120 def __init__(self, model, **kwargs): 66 121 self._sasmodel = model 67 122 pars, pd_types = create_parameters(model.info, **kwargs) 68 for k, v in pars.items():123 for k, v in pars.items(): 69 124 setattr(self, k, v) 70 for k, v in pd_types.items():125 for k, v in pd_types.items(): 71 126 setattr(self, k, v) 72 127 self._parameter_names = list(pars.keys()) … … 75 130 def parameters(self): 76 131 """ 77 Return a dictionary of parameters 132 Return a dictionary of parameters objects for the parameters, 133 excluding polydispersity distribution type. 78 134 """ 79 135 return dict((k, getattr(self, k)) for k in self._parameter_names) 80 136 81 137 def state(self): 138 """ 139 Return a dictionary of current values for all the parameters, 140 including polydispersity distribution type. 141 """ 82 142 pars = dict((k, getattr(self, k).value) for k in self._parameter_names) 83 143 pars.update((k, getattr(self, k)) for k in self._pd_type_names) … … 85 145 86 146 class Experiment(DataMixin): 87 """ 88 Return a bumps wrapper for a SAS model. 89 90 *data* is the data to be fitted. 91 92 *model* is the SAS model from :func:`core.load_model`. 147 r""" 148 Bumps wrapper for a SAS experiment. 149 150 *data* is a :class:`data.Data1D`, :class:`data.Data2D` or 151 :class:`data.Sesans` object. Use :func:`data.empty_data1D` or 152 :func:`data.empty_data2D` to define $q, \Delta q$ calculation 153 points for displaying the SANS curve when there is no measured data. 154 155 *model* is a :class:`Model` object. 93 156 94 157 *cutoff* is the integration cutoff, which avoids computing the 95 158 the SAS model where the polydispersity weight is low. 96 159 97 Model parameters can be initialized with additional keyword 98 arguments, or by assigning to model.parameter_name.value. 99 100 The resulting bumps model can be used directly in a FitProblem call. 160 The resulting model can be used directly in a Bumps FitProblem call. 101 161 """ 102 162 def __init__(self, data, model, cutoff=1e-5): … … 109 169 110 170 def update(self): 171 """ 172 Call when model parameters have changed and theory needs to be 173 recalculated. 174 """ 111 175 self._cache = {} 112 176 113 177 def numpoints(self): 114 178 """ 115 Return the number ofpoints179 Return the number of data points 116 180 """ 117 181 return len(self.Iq) … … 124 188 125 189 def theory(self): 190 """ 191 Return the theory corresponding to the model parameters. 192 193 This method uses lazy evaluation, and requires model.update() to be 194 called when the parameters have changed. 195 """ 126 196 if 'theory' not in self._cache: 127 197 pars = self.model.state() … … 130 200 131 201 def residuals(self): 202 """ 203 Return theory minus data normalized by uncertainty. 204 """ 132 205 #if np.any(self.err ==0): print("zeros in err") 133 206 return (self.theory() - self.Iq) / self.dIq 134 207 135 208 def nllf(self): 209 """ 210 Return the negative log likelihood of seeing data given the model 211 parameters, up to a normalizing constant which depends on the data 212 uncertainty. 213 """ 136 214 delta = self.residuals() 137 215 #if np.any(np.isnan(R)): print("NaN in residuals") … … 149 227 150 228 def simulate_data(self, noise=None): 229 """ 230 Generate simulated data. 231 """ 151 232 Iq = self.theory() 152 233 self._set_data(Iq, noise) 153 234 154 235 def save(self, basename): 236 """ 237 Save the model parameters and data into a file. 238 239 Not Implemented. 240 """ 155 241 pass 156 242 -
sasmodels/compare.py
r608e31e rd5e650d 28 28 29 29 from __future__ import print_function 30 31 import sys 32 import math 33 from os.path import basename, dirname, join as joinpath 34 import glob 35 import datetime 36 import traceback 37 38 import numpy as np 39 40 from . import core 41 from . import kerneldll 42 from . import generate 43 from .data import plot_theory, empty_data1D, empty_data2D 44 from .direct_model import DirectModel 45 from .convert import revert_model, constrain_new_to_old 30 46 31 47 USAGE = """ … … 72 88 # doc string so that we can display it at run time if there is an error. 73 89 # lin 74 __doc__ = __doc__ + """ 90 __doc__ = (__doc__ # pylint: disable=redefined-builtin 91 + """ 75 92 Program description 76 93 ------------------- 77 94 78 """ + USAGE 79 80 81 82 import sys 83 import math 84 from os.path import basename, dirname, join as joinpath 85 import glob 86 import datetime 87 import traceback 88 89 import numpy as np 90 95 """ 96 + USAGE) 97 98 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 99 100 # List of available models 91 101 ROOT = dirname(__file__) 92 sys.path.insert(0, ROOT) # Make sure sasmodels is first on the path93 94 95 from . import core96 from . import kerneldll97 from . import generate98 from .data import plot_theory, empty_data1D, empty_data2D99 from .direct_model import DirectModel100 from .convert import revert_model, constrain_new_to_old101 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True102 103 # List of available models104 102 MODELS = [basename(f)[:-3] 105 103 for f in sorted(glob.glob(joinpath(ROOT, "models", "[a-zA-Z]*.py")))] … … 115 113 return dt.total_seconds() 116 114 115 116 class push_seed(object): 117 """ 118 Set the seed value for the random number generator. 119 120 When used in a with statement, the random number generator state is 121 restored after the with statement is complete. 122 123 :Parameters: 124 125 *seed* : int or array_like, optional 126 Seed for RandomState 127 128 :Example: 129 130 Seed can be used directly to set the seed:: 131 132 >>> from numpy.random import randint 133 >>> push_seed(24) 134 <...push_seed object at...> 135 >>> print(randint(0,1000000,3)) 136 [242082 899 211136] 137 138 Seed can also be used in a with statement, which sets the random 139 number generator state for the enclosed computations and restores 140 it to the previous state on completion:: 141 142 >>> with push_seed(24): 143 ... print(randint(0,1000000,3)) 144 [242082 899 211136] 145 146 Using nested contexts, we can demonstrate that state is indeed 147 restored after the block completes:: 148 149 >>> with push_seed(24): 150 ... print(randint(0,1000000)) 151 ... with push_seed(24): 152 ... print(randint(0,1000000,3)) 153 ... print(randint(0,1000000)) 154 242082 155 [242082 899 211136] 156 899 157 158 The restore step is protected against exceptions in the block:: 159 160 >>> with push_seed(24): 161 ... print(randint(0,1000000)) 162 ... try: 163 ... with push_seed(24): 164 ... print(randint(0,1000000,3)) 165 ... raise Exception() 166 ... except: 167 ... print("Exception raised") 168 ... print(randint(0,1000000)) 169 242082 170 [242082 899 211136] 171 Exception raised 172 899 173 """ 174 def __init__(self, seed=None): 175 self._state = np.random.get_state() 176 np.random.seed(seed) 177 178 def __enter__(self): 179 return None 180 181 def __exit__(self, *args): 182 np.random.set_state(self._state) 117 183 118 184 def tic(): … … 177 243 return [0, (2*v if v > 0 else 1)] 178 244 245 179 246 def _randomize_one(p, v): 180 247 """ … … 186 253 return np.random.uniform(*parameter_range(p, v)) 187 254 255 188 256 def randomize_pars(pars, seed=None): 189 257 """ … … 195 263 :func:`constrain_pars` needs to be called afterward.. 196 264 """ 197 np.random.seed(seed)198 # Note: the sort guarantees order `of calls to random number generator199 pars = dict((p, _randomize_one(p, v))200 for p, v in sorted(pars.items()))265 with push_seed(seed): 266 # Note: the sort guarantees order `of calls to random number generator 267 pars = dict((p, _randomize_one(p, v)) 268 for p, v in sorted(pars.items())) 201 269 return pars 202 270 … … 281 349 theory = lambda: smearer.get_value() 282 350 else: 283 theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]]) 351 theory = lambda: model.evalDistribution([data.qx_data[index], 352 data.qy_data[index]]) 284 353 elif smearer is not None: 285 354 theory = lambda: smearer(model.evalDistribution(data.x)) … … 416 485 try: 417 486 base_value, base_time = time_calculation(base, pars, Nbase) 418 print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value))) 487 print("%s t=%.1f ms, intensity=%.0f" 488 % (base.engine, base_time, sum(base_value))) 419 489 except ImportError: 420 490 traceback.print_exc() … … 426 496 try: 427 497 comp_value, comp_time = time_calculation(comp, pars, Ncomp) 428 print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value))) 498 print("%s t=%.1f ms, intensity=%.0f" 499 % (comp.engine, comp_time, sum(comp_value))) 429 500 except ImportError: 430 501 traceback.print_exc() … … 433 504 # Compare, but only if computing both forms 434 505 if Nbase > 0 and Ncomp > 0: 435 #print("speedup %.2g"%(comp_time/base_time))436 #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value)))437 #comp *= max(base_value/comp_value)438 506 resid = (base_value - comp_value) 439 507 relerr = resid/comp_value 440 _print_stats("|%s-%s|"%(base.engine, comp.engine) + (" "*(3+len(comp.engine))), 508 _print_stats("|%s-%s|" 509 % (base.engine, comp.engine) + (" "*(3+len(comp.engine))), 441 510 resid) 442 _print_stats("|(%s-%s)/%s|"%(base.engine, comp.engine, comp.engine), 511 _print_stats("|(%s-%s)/%s|" 512 % (base.engine, comp.engine, comp.engine), 443 513 relerr) 444 514 … … 459 529 if Nbase > 0: 460 530 if Ncomp > 0: plt.subplot(131) 461 plot_theory(data, base_value, view=view, plot_data=False, limits=limits)531 plot_theory(data, base_value, view=view, use_data=False, limits=limits) 462 532 plt.title("%s t=%.1f ms"%(base.engine, base_time)) 463 533 #cbar_title = "log I" 464 534 if Ncomp > 0: 465 535 if Nbase > 0: plt.subplot(132) 466 plot_theory(data, comp_value, view=view, plot_data=False, limits=limits)536 plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 467 537 plt.title("%s t=%.1f ms"%(comp.engine, comp_time)) 468 538 #cbar_title = "log I" 469 539 if Ncomp > 0 and Nbase > 0: 470 540 plt.subplot(133) 471 if '-abs' in opts:541 if not opts['rel_err']: 472 542 err, errstr, errview = resid, "abs err", "linear" 473 543 else: 474 544 err, errstr, errview = abs(relerr), "rel err", "log" 475 545 #err,errstr = base/comp,"ratio" 476 plot_theory(data, None, resid=err, view=errview, plot_data=False) 546 plot_theory(data, None, resid=err, view=errview, use_data=False) 547 if view == 'linear': 548 plt.xscale('linear') 477 549 plt.title("max %s = %.3g"%(errstr, max(abs(err)))) 478 550 #cbar_title = errstr if errview=="linear" else "log "+errstr … … 534 606 def columnize(L, indent="", width=79): 535 607 """ 536 Format a list of strings into columns for printing. 608 Format a list of strings into columns. 609 610 Returns a string with carriage returns ready for printing. 537 611 """ 538 612 column_width = max(len(w) for w in L) + 1 … … 591 665 invalid = [o[1:] for o in flags 592 666 if o[1:] not in NAME_OPTIONS 593 667 and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 594 668 if invalid: 595 669 print("Invalid options: %s"%(", ".join(invalid))) … … 597 671 598 672 673 # pylint: disable=bad-whitespace 599 674 # Interpret the flags 600 675 opts = { … … 651 726 elif arg == '-sasview': engines.append(arg[1:]) 652 727 elif arg == '-edit': opts['explore'] = True 728 # pylint: enable=bad-whitespace 653 729 654 730 if len(engines) == 0: … … 675 751 presets = {} 676 752 for arg in values: 677 k, v = arg.split('=',1)753 k, v = arg.split('=', 1) 678 754 if k not in pars: 679 755 # extract base name without polydispersity info 680 756 s = set(p.split('_pd')[0] for p in pars) 681 print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))757 print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 682 758 sys.exit(1) 683 759 presets[k] = float(v) if not k.endswith('type') else v … … 697 773 698 774 # Create the computational engines 699 data, _ index= make_data(opts)775 data, _ = make_data(opts) 700 776 if n1: 701 777 base = make_engine(model_definition, data, engines[0], opts['cutoff']) … … 707 783 comp = None 708 784 785 # pylint: disable=bad-whitespace 709 786 # Remember it all 710 787 opts.update({ … … 718 795 'engines' : [base, comp], 719 796 }) 797 # pylint: enable=bad-whitespace 720 798 721 799 return opts 722 800 723 def main():724 opts = parse_opts()725 if opts['explore']:726 explore(opts)727 else:728 compare(opts)729 730 801 def explore(opts): 802 """ 803 Explore the model using the Bumps GUI. 804 """ 731 805 import wx 732 806 from bumps.names import FitProblem … … 734 808 735 809 problem = FitProblem(Explore(opts)) 736 is Mac = "cocoa" in wx.version()810 is_mac = "cocoa" in wx.version() 737 811 app = wx.App() 738 812 frame = AppFrame(parent=None, title="explore") 739 if not is Mac: frame.Show()813 if not is_mac: frame.Show() 740 814 frame.panel.set_model(model=problem) 741 815 frame.panel.Layout() 742 816 frame.panel.aui.Split(0, wx.TOP) 743 if is Mac: frame.Show()817 if is_mac: frame.Show() 744 818 app.MainLoop() 745 819 746 820 class Explore(object): 747 821 """ 748 Return a bumps wrapper for a SAS model comparison. 822 Bumps wrapper for a SAS model comparison. 823 824 The resulting object can be used as a Bumps fit problem so that 825 parameters can be adjusted in the GUI, with plots updated on the fly. 749 826 """ 750 827 def __init__(self, opts): … … 787 864 Return cost. 788 865 """ 866 # pylint: disable=no-self-use 789 867 return 0. # No nllf 790 868 … … 804 882 805 883 884 def main(): 885 """ 886 Main program. 887 """ 888 opts = parse_opts() 889 if opts['explore']: 890 explore(opts) 891 else: 892 compare(opts) 893 806 894 if __name__ == "__main__": 807 895 main() -
sasmodels/compare_many.py
r6458608 r4f2478e 1 1 #!/usr/bin/env python 2 """ 3 Program to compare results from many random parameter sets for a given model. 4 5 The result is a comma separated value (CSV) table that can be redirected 6 from standard output into a file and loaded into a spreadsheet. 7 8 The models are compared for each parameter set and if the difference is 9 greater than expected for that precision, the parameter set is labeled 10 as bad and written to the output, along with the random seed used to 11 generate that parameter value. This seed can be used with :mod:`compare` 12 to reload and display the details of the model. 13 """ 2 14 from __future__ import print_function 3 15 … … 14 26 15 27 def calc_stats(target, value, index): 28 """ 29 Calculate statistics between the target value and the computed value. 30 31 *target* and *value* are the vectors being compared, with the 32 difference normalized by *target* to get relative error. Only 33 the elements listed in *index* are used, though index may be 34 and empty slice defined by *slice(None, None)*. 35 36 Returns: 37 38 *maxrel* the maximum relative difference 39 40 *rel95* the relative difference with the 5% biggest differences ignored 41 42 *maxabs* the maximum absolute difference for the 5% biggest differences 43 44 *maxval* the maximum value for the 5% biggest differences 45 """ 16 46 resid = abs(value-target)[index] 17 47 relerr = resid/target[index] 18 s rel= np.argsort(relerr)48 sorted_rel_index = np.argsort(relerr) 19 49 #p90 = int(len(relerr)*0.90) 20 50 p95 = int(len(relerr)*0.95) 21 51 maxrel = np.max(relerr) 22 rel95 = relerr[s rel[p95]]23 maxabs = np.max(resid[s rel[p95:]])24 maxval = np.max(value[s rel[p95:]])25 return maxrel, rel95,maxabs,maxval52 rel95 = relerr[sorted_rel_index[p95]] 53 maxabs = np.max(resid[sorted_rel_index[p95:]]) 54 maxval = np.max(value[sorted_rel_index[p95:]]) 55 return maxrel, rel95, maxabs, maxval 26 56 27 57 def print_column_headers(pars, parts): 58 """ 59 Generate column headers for the differences and for the parameters, 60 and print them to standard output. 61 """ 28 62 stats = list('Max rel err|95% rel err|Max abs err above 90% rel|Max value above 90% rel'.split('|')) 29 63 groups = [''] … … 36 70 print(','.join('"%s"'%c for c in columns)) 37 71 72 # Target 'good' value for various precision levels. 38 73 PRECISION = { 39 74 'fast': 1e-3, … … 48 83 def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5, 49 84 base='sasview', comp='double'): 50 is2D = hasattr(data, 'qx_data') 85 r""" 86 Compare the model under different calculation engines. 87 88 *name* is the name of the model. 89 90 *data* is the data object giving $q, \Delta q$ calculation points. 91 92 *index* is the active set of points. 93 94 *N* is the number of comparisons to make. 95 96 *cutoff* is the polydispersity weight cutoff to make the calculation 97 a little bit faster. 98 99 *base* and *comp* are the names of the calculation engines to compare. 100 """ 101 102 is_2d = hasattr(data, 'qx_data') 51 103 model_definition = core.load_model_definition(name) 52 104 pars = get_demo_pars(model_definition) 53 105 header = ('\n"Model","%s","Count","%d","Dimension","%s"' 54 % (name, N, "2D" if is 2Delse "1D"))106 % (name, N, "2D" if is_2d else "1D")) 55 107 if not mono: header += ',"Cutoff",%g'%(cutoff,) 56 108 print(header) 57 109 58 if is 2D:110 if is_2d: 59 111 info = generate.make_info(model_definition) 60 112 partype = info['partype'] … … 69 121 # declarations are not available in python 2.7. 70 122 def try_model(fn, pars): 123 """ 124 Return the model evaluated at *pars*. If there is an exception, 125 print it and return NaN of the right shape. 126 """ 71 127 try: 72 128 result = fn(**pars) … … 82 138 return result 83 139 def check_model(pars): 140 """ 141 Run the two calculators against *pars*, returning statistics 142 on the differences. See :func:`calc_stats` for the list of stats. 143 """ 84 144 base_value = try_model(calc_base, pars) 85 145 comp_value = try_model(calc_comp, pars) … … 108 168 good = [True] 109 169 columns = check_model(pars_i) 110 columns += [v for _, v in sorted(pars_i.items())]170 columns += [v for _, v in sorted(pars_i.items())] 111 171 if first: 112 172 labels = [" vs. ".join((calc_base.engine, calc_comp.engine))] … … 121 181 122 182 def print_usage(): 183 """ 184 Print the command usage string. 185 """ 123 186 print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)") 124 187 125 188 126 189 def print_models(): 190 """ 191 Print the list of available models in columns. 192 """ 127 193 print(columnize(MODELS, indent=" ")) 128 194 129 195 130 196 def print_help(): 197 """ 198 Print usage string, the option description and the list of available models. 199 """ 131 200 print_usage() 132 201 print("""\ … … 158 227 159 228 def main(): 160 if len(sys.argv) not in (6,7): 229 """ 230 Main program. 231 """ 232 if len(sys.argv) not in (6, 7): 161 233 print_help() 162 234 sys.exit(1) … … 182 254 183 255 data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 184 256 'accuracy': 'Low', 'view':'log'}) 185 257 model_list = [model] if model != "all" else MODELS 186 258 for model in model_list: … … 189 261 190 262 if __name__ == "__main__": 263 #from .compare import push_seed 264 #with push_seed(1): main() 191 265 main() -
sasmodels/convert.py
r03582f9 rd368d21 1 """ 2 Convert models to and from sasview. 3 """ 4 import warnings 5 6 PD_DOT = [ 7 ("", ""), 8 ("_pd", ".width"), 9 ("_pd_n", ".npts"), 10 ("_pd_nsigma", ".nsigmas"), 11 ("_pd_type", ".type"), 12 ] 13 def _convert_pars(pars, mapping): 14 """ 15 Rename the parameters and any associated polydispersity attributes. 16 """ 17 newpars = pars.copy() 18 for new, old in mapping.items(): 19 if old == new: continue 20 for pd, dot in PD_DOT: 21 if old+dot in newpars: 22 if new is not None: 23 newpars[new+pd] = pars[old+dot] 24 del newpars[old+dot] 25 return newpars 26 27 def _rescale_sld(pars): 28 """ 29 rescale all sld parameters in the new model definition by 1e6 so the 30 numbers are nicer. Relies on the fact that all sld parameters in the 31 new model definition end with sld. 32 """ 33 return dict((p, (v*1e6 if p.endswith('sld') else v)) 34 for p,v in pars.items()) 35 36 def convert_model(name, pars): 37 """ 38 Convert model from old style parameter names to new style. 39 """ 40 _,_ = name,pars # lint 41 raise NotImplementedError 42 # need to load all new models in order to determine old=>new 43 # model name mapping 44 45 def _unscale_sld(pars): 46 """ 47 rescale all sld parameters in the new model definition by 1e6 so the 48 numbers are nicer. Relies on the fact that all sld parameters in the 49 new model definition end with sld. 50 """ 51 return dict((p, (v*1e-6 if p.endswith('sld') else v)) 52 for p,v in pars.items()) 53 54 def _remove_pd(pars, key, name): 55 """ 56 Remove polydispersity from the parameter list. 57 58 Note: operates in place 59 """ 60 # Bumps style parameter names 61 pd = pars.pop(key+".width", 0.0) 62 pd_n = pars.pop(key+".npts", 0) 63 if pd != 0.0 and pd_n != 0: 64 warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name)) 65 pars.pop(key+".nsigmas", None) 66 pars.pop(key+".type", None) 67 return pars 68 69 def _revert_pars(pars, mapping): 70 """ 71 Rename the parameters and any associated polydispersity attributes. 72 """ 73 newpars = pars.copy() 74 75 for new, old in mapping.items(): 76 for pd, dot in PD_DOT: 77 if old and old+pd == new+dot: 78 continue 79 if new+pd in newpars: 80 if old is not None: 81 newpars[old+dot] = pars[new+pd] 82 del newpars[new+pd] 83 for k in list(newpars.keys()): 84 for pd, dot in PD_DOT[1:]: # skip "" => "" 85 if k.endswith(pd): 86 newpars[k[:-len(pd)]+dot] = newpars[k] 87 del newpars[k] 88 return newpars 89 90 def revert_model(model_definition, pars): 91 """ 92 Convert model from new style parameter names to old style. 93 """ 94 mapping = model_definition.oldpars 95 oldname = model_definition.oldname 96 oldpars = _revert_pars(_unscale_sld(pars), mapping) 97 98 # Note: update compare.constrain_pars to match 99 name = model_definition.name 100 if name in ('teubner_strey', 'broad_peak', 'two_lorentzian', 'gel_fit', 101 'binary_hard_sphere'): 102 if oldpars.pop('scale', 1.0) != 1.0: 103 warnings.warn("parameter scale not used in sasview %s"%name) 104 elif name in ('guinier',): 105 if oldpars.pop('background', 0.0) != 0.0: 106 warnings.warn("parameter background not used in sasview %s"%name) 107 elif getattr(model_definition, 'category', None) == 'structure-factor': 108 if oldpars.pop('scale', 1.0) != 1.0: 109 warnings.warn("parameter scale not used in sasview %s"%name) 110 if oldpars.pop('background', 0.0) != 0.0: 111 warnings.warn("parameter background not used in sasview %s"%name) 112 elif name == 'pearl_necklace': 113 _remove_pd(oldpars, 'num_pearls', name) 114 _remove_pd(oldpars, 'thick_string', name) 115 elif name == 'rpa': 116 # convert scattering lengths from femtometers to centimeters 117 for p in "La", "Lb", "Lc", "Ld": 118 if p in oldpars: oldpars[p] *= 1e-13 119 120 return oldname, oldpars 121 122 def constrain_new_to_old(model_definition, pars): 123 """ 124 Restrict parameter values to those that will match sasview. 125 """ 126 # Note: update convert.revert_model to match 127 name = model_definition.name 128 if name in ('teubner_strey','broad_peak', 'two_lorentzian', 'gel_fit'): 129 pars['scale'] = 1 130 elif name in ('guinier',): 131 pars['background'] = 0 132 elif name == 'pearl_necklace': 133 pars['string_thickness_pd_n'] = 0 134 pars['number_of_pearls_pd_n'] = 0 135 elif name == 'rpa': 136 pars['case_num'] = int(pars['case_num']) 137 elif getattr(model_definition, 'category', None) == 'structure-factor': 138 pars['scale'], pars['background'] = 1, 0 139 1 2 3 4 5 <!DOCTYPE html> 6 <html lang="en" class=" is-copy-enabled is-u2f-enabled"> 7 <head prefix="og: http://ogp.me/ns# fb: http://ogp.me/ns/fb# object: http://ogp.me/ns/object# article: http://ogp.me/ns/article# profile: http://ogp.me/ns/profile#"> 8 <meta charset='utf-8'> 9 10 <link crossorigin="anonymous" href="https://assets-cdn.github.com/assets/github-e3dd2ae433414e240167f740f90ff2599e28804648787933de7c0183fae81c29.css" integrity="sha256-490q5DNBTiQBZ/dA+Q/yWZ4ogEZIeHkz3nwBg/roHCk=" media="all" rel="stylesheet" /> 11 <link crossorigin="anonymous" href="https://assets-cdn.github.com/assets/github2-fd8d48abb9063f51f186b0da98caa88d32b7ca8baecaa10d8a91c18a6f129b7f.css" integrity="sha256-/Y1Iq7kGP1HxhrDamMqojTK3youuyqENipHBim8Sm38=" media="all" rel="stylesheet" /> 12 13 14 15 16 <link as="script" href="https://assets-cdn.github.com/assets/frameworks-ee521b8e9facac68ff27e93fc3ae0f8ed811d7bf9e434e84f4b9ea227780b084.js" rel="preload" /> 17 <link as="script" href="https://assets-cdn.github.com/assets/github-696336964b7c42e4c6f4dfbaf1f8e57f425cd9a9d18a12f3fe6e3dd744fd7d13.js" rel="preload" /> 18 19 <meta http-equiv="X-UA-Compatible" content="IE=edge"> 20 <meta http-equiv="Content-Language" content="en"> 21 <meta name="viewport" content="width=1020"> 22 23 24 <title>sasmodels/convert.py at master · SasView/sasmodels · GitHub</title> 25 <link rel="search" type="application/opensearchdescription+xml" href="/opensearch.xml" title="GitHub"> 26 <link rel="fluid-icon" href="https://github.com/fluidicon.png" title="GitHub"> 27 <link rel="apple-touch-icon" href="/apple-touch-icon.png"> 28 <link rel="apple-touch-icon" sizes="57x57" href="/apple-touch-icon-57x57.png"> 29 <link rel="apple-touch-icon" sizes="60x60" href="/apple-touch-icon-60x60.png"> 30 <link rel="apple-touch-icon" sizes="72x72" href="/apple-touch-icon-72x72.png"> 31 <link rel="apple-touch-icon" sizes="76x76" href="/apple-touch-icon-76x76.png"> 32 <link rel="apple-touch-icon" sizes="114x114" href="/apple-touch-icon-114x114.png"> 33 <link rel="apple-touch-icon" sizes="120x120" href="/apple-touch-icon-120x120.png"> 34 <link rel="apple-touch-icon" sizes="144x144" href="/apple-touch-icon-144x144.png"> 35 <link rel="apple-touch-icon" sizes="152x152" href="/apple-touch-icon-152x152.png"> 36 <link rel="apple-touch-icon" sizes="180x180" href="/apple-touch-icon-180x180.png"> 37 <meta property="fb:app_id" content="1401488693436528"> 38 39 <meta content="https://avatars1.githubusercontent.com/u/10882470?v=3&s=400" name="twitter:image:src" /><meta content="@github" name="twitter:site" /><meta content="summary" name="twitter:card" /><meta content="SasView/sasmodels" name="twitter:title" /><meta content="sasmodels - Package for calculation of small angle scattering models using OpenCL. Builds here: https://jenkins.esss.dk/sasview/view/Sasmodels-Builds/" name="twitter:description" /> 40 <meta content="https://avatars1.githubusercontent.com/u/10882470?v=3&s=400" property="og:image" /><meta content="GitHub" property="og:site_name" /><meta content="object" property="og:type" /><meta content="SasView/sasmodels" property="og:title" /><meta content="https://github.com/SasView/sasmodels" property="og:url" /><meta content="sasmodels - Package for calculation of small angle scattering models using OpenCL. Builds here: https://jenkins.esss.dk/sasview/view/Sasmodels-Builds/" property="og:description" /> 41 <meta name="browser-stats-url" content="https://api.github.com/_private/browser/stats"> 42 <meta name="browser-errors-url" content="https://api.github.com/_private/browser/errors"> 43 <link rel="assets" href="https://assets-cdn.github.com/"> 44 45 <meta name="pjax-timeout" content="1000"> 46 47 48 <meta name="msapplication-TileImage" content="/windows-tile.png"> 49 <meta name="msapplication-TileColor" content="#ffffff"> 50 <meta name="selected-link" value="repo_source" data-pjax-transient> 51 52 <meta name="google-site-verification" content="KT5gs8h0wvaagLKAVWq8bbeNwnZZK1r1XQysX3xurLU"> 53 <meta name="google-site-verification" content="ZzhVyEFwb7w3e0-uOTltm8Jsck2F5StVihD0exw2fsA"> 54 <meta name="google-analytics" content="UA-3769691-2"> 55 56 <meta content="collector.githubapp.com" name="octolytics-host" /><meta content="github" name="octolytics-app-id" /><meta content="32FD15EB:517B:95519E0:56B61366" name="octolytics-dimension-request_id" /> 57 <meta content="/<user-name>/<repo-name>/blob/show" data-pjax-transient="true" name="analytics-location" /> 58 59 60 61 <meta class="js-ga-set" name="dimension1" content="Logged Out"> 62 63 64 65 <meta name="hostname" content="github.com"> 66 <meta name="user-login" content=""> 67 68 <meta name="expected-hostname" content="github.com"> 69 70 <link rel="mask-icon" href="https://assets-cdn.github.com/pinned-octocat.svg" color="#4078c0"> 71 <link rel="icon" type="image/x-icon" href="https://assets-cdn.github.com/favicon.ico"> 72 73 <meta content="ed1699ee8f88e1fc0f2738f860d9b53cf2d34560" name="form-nonce" /> 74 75 <meta http-equiv="x-pjax-version" content="2330245d99bfcba85ea2b6d282c0cb53"> 76 77 78 <meta name="description" content="sasmodels - Package for calculation of small angle scattering models using OpenCL. Builds here: https://jenkins.esss.dk/sasview/view/Sasmodels-Builds/"> 79 <meta name="go-import" content="github.com/SasView/sasmodels git https://github.com/SasView/sasmodels.git"> 80 81 <meta content="10882470" name="octolytics-dimension-user_id" /><meta content="SasView" name="octolytics-dimension-user_login" /><meta content="30761174" name="octolytics-dimension-repository_id" /><meta content="SasView/sasmodels" name="octolytics-dimension-repository_nwo" /><meta content="true" name="octolytics-dimension-repository_public" /><meta content="false" name="octolytics-dimension-repository_is_fork" /><meta content="30761174" name="octolytics-dimension-repository_network_root_id" /><meta content="SasView/sasmodels" name="octolytics-dimension-repository_network_root_nwo" /> 82 <link href="https://github.com/SasView/sasmodels/commits/master.atom" rel="alternate" title="Recent Commits to sasmodels:master" type="application/atom+xml"> 83 84 85 <link rel="canonical" href="https://github.com/SasView/sasmodels/blob/master/sasmodels/convert.py" data-pjax-transient> 86 </head> 87 88 89 <body class="logged_out env-production windows vis-public page-blob"> 90 <a href="#start-of-content" tabindex="1" class="accessibility-aid js-skip-to-content">Skip to content</a> 91 92 93 94 95 96 97 98 99 <div class="header header-logged-out" role="banner"> 100 <div class="container clearfix"> 101 102 <a class="header-logo-wordmark" href="https://github.com/" data-ga-click="(Logged out) Header, go to homepage, icon:logo-wordmark"> 103 <svg aria-hidden="true" class="octicon octicon-logo-github" height="28" role="img" version="1.1" viewBox="0 0 45 16" width="78"><path d="M8.64 5.19H4.88c-0.11 0-0.19 0.08-0.19 0.17v1.84c0 0.09 0.08 0.17 0.19 0.17h1.47v2.3s-0.33 0.11-1.25 0.11c-1.08 0-2.58-0.39-2.58-3.7s1.58-3.73 3.05-3.73c1.27 0 1.81 0.22 2.17 0.33 0.11 0.03 0.2-0.08 0.2-0.17l0.42-1.78c0-0.05-0.02-0.09-0.06-0.14-0.14-0.09-1.02-0.58-3.2-0.58C2.58 0 0 1.06 0 6.2s2.95 5.92 5.44 5.92c2.06 0 3.31-0.89 3.31-0.89 0.05-0.02 0.06-0.09 0.06-0.13V5.36c0-0.09-0.08-0.17-0.19-0.17h0.02zM27.7 0.44h-2.13c-0.09 0-0.17 0.08-0.17 0.17v4.09h-3.31V0.61c0-0.09-0.08-0.17-0.17-0.17h-2.13c-0.09 0-0.17 0.08-0.17 0.17v11.11c0 0.09 0.09 0.17 0.17 0.17h2.13c0.09 0 0.17-0.08 0.17-0.17V6.97h3.31l-0.02 4.75c0 0.09 0.08 0.17 0.17 0.17h2.13c0.09 0 0.17-0.08 0.17-0.17V0.61c0-0.09-0.08-0.17-0.17-0.17h0.02zM11.19 0.69c-0.77 0-1.38 0.61-1.38 1.38s0.61 1.38 1.38 1.38c0.75 0 1.36-0.61 1.36-1.38s-0.61-1.38-1.36-1.38z m1.22 3.55c0-0.09-0.08-0.17-0.17-0.17H10.11c-0.09 0-0.17 0.09-0.17 0.2 0 0 0 6.17 0 7.34 0 0.2 0.13 0.27 0.3 0.27 0 0 0.91 0 1.92 0 0.2 0 0.25-0.09 0.25-0.27 0-0.39 0-7.36 0-7.36v-0.02z m23.52-0.16h-2.09c-0.11 0-0.17 0.08-0.17 0.19v5.44s-0.55 0.39-1.3 0.39-0.97-0.34-0.97-1.09c0-0.73 0-4.75 0-4.75 0-0.09-0.08-0.17-0.17-0.17h-2.14c-0.09 0-0.17 0.08-0.17 0.17 0 0 0 2.91 0 5.11s1.23 2.75 2.92 2.75c1.39 0 2.52-0.77 2.52-0.77s0.05 0.39 0.08 0.45c0.02 0.05 0.09 0.09 0.16 0.09h1.34c0.11 0 0.17-0.08 0.17-0.17l0.02-7.47c0-0.09-0.08-0.17-0.19-0.17z m5.77-0.25c-1.2 0-2.02 0.53-2.02 0.53V0.59c0-0.09-0.08-0.17-0.17-0.17h-2.13c-0.09 0-0.17 0.08-0.17 0.17l-0.02 11.11c0 0.09 0.09 0.17 0.19 0.17h1.48c0.06 0 0.11-0.02 0.14-0.08 0.05-0.06 0.09-0.52 0.09-0.52s0.88 0.83 2.52 0.83c1.94 0 3.05-0.98 3.05-4.41s-1.77-3.88-2.97-3.88z m-0.83 6.27c-0.73-0.02-1.22-0.36-1.22-0.36V6.22s0.48-0.3 1.08-0.34c0.77-0.08 1.5 0.16 1.5 1.97 0 1.91-0.33 2.28-1.36 2.25z m-22.33-0.05c-0.09 0-0.33 0.05-0.58 0.05-0.78 0-1.05-0.36-1.05-0.83s0-3.13 0-3.13h1.59c0.09 0 0.16-0.08 0.16-0.19V4.25c0-0.09-0.08-0.17-0.16-0.17h-1.59V1.97c0-0.08-0.05-0.13-0.14-0.13H14.61c-0.09 0-0.14 0.05-0.14 0.13v2.17s-1.09 0.27-1.16 0.28c-0.08 0.02-0.13 0.09-0.13 0.17v1.36c0 0.11 0.08 0.19 0.17 0.19h1.11s0 1.44 0 3.28c0 2.44 1.7 2.69 2.86 2.69 0.53 0 1.17-0.17 1.27-0.22 0.06-0.02 0.09-0.09 0.09-0.16v-1.5c0-0.11-0.08-0.19-0.17-0.19h0.02z"></path></svg> 104 </a> 105 106 <div class="header-actions" role="navigation"> 107 <a class="btn btn-primary" href="/join?source=header-repo" data-ga-click="(Logged out) Header, clicked Sign up, text:sign-up">Sign up</a> 108 <a class="btn" href="/login?return_to=%2FSasView%2Fsasmodels%2Fblob%2Fmaster%2Fsasmodels%2Fconvert.py" data-ga-click="(Logged out) Header, clicked Sign in, text:sign-in">Sign in</a> 109 </div> 110 111 <div class="site-search repo-scope js-site-search" role="search"> 112 <!-- </textarea> --><!-- '"` --><form accept-charset="UTF-8" action="/SasView/sasmodels/search" class="js-site-search-form" data-global-search-url="/search" data-repo-search-url="/SasView/sasmodels/search" method="get"><div style="margin:0;padding:0;display:inline"><input name="utf8" type="hidden" value="✓" /></div> 113 <label class="js-chromeless-input-container form-control"> 114 <div class="scope-badge">This repository</div> 115 <input type="text" 116 class="js-site-search-focus js-site-search-field is-clearable chromeless-input" 117 data-hotkey="s" 118 name="q" 119 placeholder="Search" 120 aria-label="Search this repository" 121 data-global-scope-placeholder="Search GitHub" 122 data-repo-scope-placeholder="Search" 123 tabindex="1" 124 autocapitalize="off"> 125 </label> 126 </form> 127 </div> 128 129 <ul class="header-nav left" role="navigation"> 130 <li class="header-nav-item"> 131 <a class="header-nav-link" href="/explore" data-ga-click="(Logged out) Header, go to explore, text:explore">Explore</a> 132 </li> 133 <li class="header-nav-item"> 134 <a class="header-nav-link" href="/features" data-ga-click="(Logged out) Header, go to features, text:features">Features</a> 135 </li> 136 <li class="header-nav-item"> 137 <a class="header-nav-link" href="https://enterprise.github.com/" data-ga-click="(Logged out) Header, go to enterprise, text:enterprise">Enterprise</a> 138 </li> 139 <li class="header-nav-item"> 140 <a class="header-nav-link" href="/pricing" data-ga-click="(Logged out) Header, go to pricing, text:pricing">Pricing</a> 141 </li> 142 </ul> 143 144 </div> 145 </div> 146 147 148 149 <div id="start-of-content" class="accessibility-aid"></div> 150 151 <div id="js-flash-container"> 152 </div> 153 154 155 <div role="main" class="main-content"> 156 <div itemscope itemtype="http://schema.org/WebPage"> 157 <div id="js-repo-pjax-container" class="context-loader-container js-repo-nav-next" data-pjax-container> 158 159 <div class="pagehead repohead instapaper_ignore readability-menu experiment-repo-nav"> 160 <div class="container repohead-details-container"> 161 162 163 164 <ul class="pagehead-actions"> 165 166 <li> 167 <a href="/login?return_to=%2FSasView%2Fsasmodels" 168 class="btn btn-sm btn-with-count tooltipped tooltipped-n" 169 aria-label="You must be signed in to watch a repository" rel="nofollow"> 170 <svg aria-hidden="true" class="octicon octicon-eye" height="16" role="img" version="1.1" viewBox="0 0 16 16" width="16"><path d="M8.06 2C3 2 0 8 0 8s3 6 8.06 6c4.94 0 7.94-6 7.94-6S13 2 8.06 2z m-0.06 10c-2.2 0-4-1.78-4-4 0-2.2 1.8-4 4-4 2.22 0 4 1.8 4 4 0 2.22-1.78 4-4 4z m2-4c0 1.11-0.89 2-2 2s-2-0.89-2-2 0.89-2 2-2 2 0.89 2 2z"></path></svg> 171 Watch 172 </a> 173 <a class="social-count" href="/SasView/sasmodels/watchers"> 174 16 175 </a> 176 177 </li> 178 179 <li> 180 <a href="/login?return_to=%2FSasView%2Fsasmodels" 181 class="btn btn-sm btn-with-count tooltipped tooltipped-n" 182 aria-label="You must be signed in to star a repository" rel="nofollow"> 183 <svg aria-hidden="true" class="octicon octicon-star" height="16" role="img" version="1.1" viewBox="0 0 14 16" width="14"><path d="M14 6l-4.9-0.64L7 1 4.9 5.36 0 6l3.6 3.26L2.67 14l4.33-2.33 4.33 2.33L10.4 9.26 14 6z"></path></svg> 184 Star 185 </a> 186 187 <a class="social-count js-social-count" href="/SasView/sasmodels/stargazers"> 188 0 189 </a> 190 191 </li> 192 193 <li> 194 <a href="/login?return_to=%2FSasView%2Fsasmodels" 195 class="btn btn-sm btn-with-count tooltipped tooltipped-n" 196 aria-label="You must be signed in to fork a repository" rel="nofollow"> 197 <svg aria-hidden="true" class="octicon octicon-repo-forked" height="16" role="img" version="1.1" viewBox="0 0 10 16" width="10"><path d="M8 1c-1.11 0-2 0.89-2 2 0 0.73 0.41 1.38 1 1.72v1.28L5 8 3 6v-1.28c0.59-0.34 1-0.98 1-1.72 0-1.11-0.89-2-2-2S0 1.89 0 3c0 0.73 0.41 1.38 1 1.72v1.78l3 3v1.78c-0.59 0.34-1 0.98-1 1.72 0 1.11 0.89 2 2 2s2-0.89 2-2c0-0.73-0.41-1.38-1-1.72V9.5l3-3V4.72c0.59-0.34 1-0.98 1-1.72 0-1.11-0.89-2-2-2zM2 4.2c-0.66 0-1.2-0.55-1.2-1.2s0.55-1.2 1.2-1.2 1.2 0.55 1.2 1.2-0.55 1.2-1.2 1.2z m3 10c-0.66 0-1.2-0.55-1.2-1.2s0.55-1.2 1.2-1.2 1.2 0.55 1.2 1.2-0.55 1.2-1.2 1.2z m3-10c-0.66 0-1.2-0.55-1.2-1.2s0.55-1.2 1.2-1.2 1.2 0.55 1.2 1.2-0.55 1.2-1.2 1.2z"></path></svg> 198 Fork 199 </a> 200 201 <a href="/SasView/sasmodels/network" class="social-count"> 202 0 203 </a> 204 </li> 205 </ul> 206 207 <h1 itemscope itemtype="http://data-vocabulary.org/Breadcrumb" class="entry-title public "> 208 <svg aria-hidden="true" class="octicon octicon-repo" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M4 9h-1v-1h1v1z m0-3h-1v1h1v-1z m0-2h-1v1h1v-1z m0-2h-1v1h1v-1z m8-1v12c0 0.55-0.45 1-1 1H6v2l-1.5-1.5-1.5 1.5V14H1c-0.55 0-1-0.45-1-1V1C0 0.45 0.45 0 1 0h10c0.55 0 1 0.45 1 1z m-1 10H1v2h2v-1h3v1h5V11z m0-10H2v9h9V1z"></path></svg> 209 <span class="author"><a href="/SasView" class="url fn" itemprop="url" rel="author"><span itemprop="title">SasView</span></a></span><!-- 210 --><span class="path-divider">/</span><!-- 211 --><strong><a href="/SasView/sasmodels" data-pjax="#js-repo-pjax-container">sasmodels</a></strong> 212 213 <span class="page-context-loader"> 214 <img alt="" height="16" src="https://assets-cdn.github.com/images/spinners/octocat-spinner-32.gif" width="16" /> 215 </span> 216 217 </h1> 218 219 </div> 220 <div class="container"> 221 222 <nav class="reponav js-repo-nav js-sidenav-container-pjax js-octicon-loaders" 223 role="navigation" 224 data-pjax="#js-repo-pjax-container"> 225 226 <a href="/SasView/sasmodels" aria-label="Code" aria-selected="true" class="js-selected-navigation-item selected reponav-item" data-hotkey="g c" data-selected-links="repo_source repo_downloads repo_commits repo_releases repo_tags repo_branches /SasView/sasmodels"> 227 <svg aria-hidden="true" class="octicon octicon-code" height="16" role="img" version="1.1" viewBox="0 0 14 16" width="14"><path d="M9.5 3l-1.5 1.5 3.5 3.5L8 11.5l1.5 1.5 4.5-5L9.5 3zM4.5 3L0 8l4.5 5 1.5-1.5L2.5 8l3.5-3.5L4.5 3z"></path></svg> 228 Code 229 </a> 230 231 <a href="/SasView/sasmodels/pulls" class="js-selected-navigation-item reponav-item" data-hotkey="g p" data-selected-links="repo_pulls /SasView/sasmodels/pulls"> 232 <svg aria-hidden="true" class="octicon octicon-git-pull-request" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M11 11.28c0-1.73 0-6.28 0-6.28-0.03-0.78-0.34-1.47-0.94-2.06s-1.28-0.91-2.06-0.94c0 0-1.02 0-1 0V0L4 3l3 3V4h1c0.27 0.02 0.48 0.11 0.69 0.31s0.3 0.42 0.31 0.69v6.28c-0.59 0.34-1 0.98-1 1.72 0 1.11 0.89 2 2 2s2-0.89 2-2c0-0.73-0.41-1.38-1-1.72z m-1 2.92c-0.66 0-1.2-0.55-1.2-1.2s0.55-1.2 1.2-1.2 1.2 0.55 1.2 1.2-0.55 1.2-1.2 1.2zM4 3c0-1.11-0.89-2-2-2S0 1.89 0 3c0 0.73 0.41 1.38 1 1.72 0 1.55 0 5.56 0 6.56-0.59 0.34-1 0.98-1 1.72 0 1.11 0.89 2 2 2s2-0.89 2-2c0-0.73-0.41-1.38-1-1.72V4.72c0.59-0.34 1-0.98 1-1.72z m-0.8 10c0 0.66-0.55 1.2-1.2 1.2s-1.2-0.55-1.2-1.2 0.55-1.2 1.2-1.2 1.2 0.55 1.2 1.2z m-1.2-8.8c-0.66 0-1.2-0.55-1.2-1.2s0.55-1.2 1.2-1.2 1.2 0.55 1.2 1.2-0.55 1.2-1.2 1.2z"></path></svg> 233 Pull requests 234 <span class="counter">0</span> 235 </a> 236 237 <a href="/SasView/sasmodels/pulse" class="js-selected-navigation-item reponav-item" data-selected-links="pulse /SasView/sasmodels/pulse"> 238 <svg aria-hidden="true" class="octicon octicon-pulse" height="16" role="img" version="1.1" viewBox="0 0 14 16" width="14"><path d="M11.5 8L8.8 5.4 6.6 8.5 5.5 1.6 2.38 8H0V10h3.6L4.5 8.2l0.9 5.4L9 8.5l1.6 1.5H14V8H11.5z"></path></svg> 239 Pulse 240 </a> 241 <a href="/SasView/sasmodels/graphs" class="js-selected-navigation-item reponav-item" data-selected-links="repo_graphs repo_contributors /SasView/sasmodels/graphs"> 242 <svg aria-hidden="true" class="octicon octicon-graph" height="16" role="img" version="1.1" viewBox="0 0 16 16" width="16"><path d="M16 14v1H0V0h1v14h15z m-11-1H3V8h2v5z m4 0H7V3h2v10z m4 0H11V6h2v7z"></path></svg> 243 Graphs 244 </a> 245 246 </nav> 247 248 </div> 249 </div> 250 251 <div class="container new-discussion-timeline experiment-repo-nav"> 252 <div class="repository-content"> 253 254 255 256 <a href="/SasView/sasmodels/blob/d5e650da432c58f8cd7c4d86e8ecb26828e17a21/sasmodels/convert.py" class="hidden js-permalink-shortcut" data-hotkey="y">Permalink</a> 257 258 <!-- blob contrib key: blob_contributors:v21:db8856a43e364dee3a304855aa1d4b41 --> 259 260 <div class="file-navigation js-zeroclipboard-container"> 261 262 <div class="select-menu js-menu-container js-select-menu left"> 263 <button class="btn btn-sm select-menu-button js-menu-target css-truncate" data-hotkey="w" 264 title="master" 265 type="button" aria-label="Switch branches or tags" tabindex="0" aria-haspopup="true"> 266 <i>Branch:</i> 267 <span class="js-select-button css-truncate-target">master</span> 268 </button> 269 270 <div class="select-menu-modal-holder js-menu-content js-navigation-container" data-pjax aria-hidden="true"> 271 272 <div class="select-menu-modal"> 273 <div class="select-menu-header"> 274 <svg aria-label="Close" class="octicon octicon-x js-menu-close" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M7.48 8l3.75 3.75-1.48 1.48-3.75-3.75-3.75 3.75-1.48-1.48 3.75-3.75L0.77 4.25l1.48-1.48 3.75 3.75 3.75-3.75 1.48 1.48-3.75 3.75z"></path></svg> 275 <span class="select-menu-title">Switch branches/tags</span> 276 </div> 277 278 <div class="select-menu-filters"> 279 <div class="select-menu-text-filter"> 280 <input type="text" aria-label="Filter branches/tags" id="context-commitish-filter-field" class="js-filterable-field js-navigation-enable" placeholder="Filter branches/tags"> 281 </div> 282 <div class="select-menu-tabs"> 283 <ul> 284 <li class="select-menu-tab"> 285 <a href="#" data-tab-filter="branches" data-filter-placeholder="Filter branches/tags" class="js-select-menu-tab" role="tab">Branches</a> 286 </li> 287 <li class="select-menu-tab"> 288 <a href="#" data-tab-filter="tags" data-filter-placeholder="Find a tagâŠ" class="js-select-menu-tab" role="tab">Tags</a> 289 </li> 290 </ul> 291 </div> 292 </div> 293 294 <div class="select-menu-list select-menu-tab-bucket js-select-menu-tab-bucket" data-tab-filter="branches" role="menu"> 295 296 <div data-filterable-for="context-commitish-filter-field" data-filterable-type="substring"> 297 298 299 <a class="select-menu-item js-navigation-item js-navigation-open " 300 href="/SasView/sasmodels/blob/develop_ric/sasmodels/convert.py" 301 data-name="develop_ric" 302 data-skip-pjax="true" 303 rel="nofollow"> 304 <svg aria-hidden="true" class="octicon octicon-check select-menu-item-icon" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M12 5L4 13 0 9l1.5-1.5 2.5 2.5 6.5-6.5 1.5 1.5z"></path></svg> 305 <span class="select-menu-item-text css-truncate-target" title="develop_ric"> 306 develop_ric 307 </span> 308 </a> 309 <a class="select-menu-item js-navigation-item js-navigation-open " 310 href="/SasView/sasmodels/blob/gh-pages/sasmodels/convert.py" 311 data-name="gh-pages" 312 data-skip-pjax="true" 313 rel="nofollow"> 314 <svg aria-hidden="true" class="octicon octicon-check select-menu-item-icon" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M12 5L4 13 0 9l1.5-1.5 2.5 2.5 6.5-6.5 1.5 1.5z"></path></svg> 315 <span class="select-menu-item-text css-truncate-target" title="gh-pages"> 316 gh-pages 317 </span> 318 </a> 319 <a class="select-menu-item js-navigation-item js-navigation-open selected" 320 href="/SasView/sasmodels/blob/master/sasmodels/convert.py" 321 data-name="master" 322 data-skip-pjax="true" 323 rel="nofollow"> 324 <svg aria-hidden="true" class="octicon octicon-check select-menu-item-icon" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M12 5L4 13 0 9l1.5-1.5 2.5 2.5 6.5-6.5 1.5 1.5z"></path></svg> 325 <span class="select-menu-item-text css-truncate-target" title="master"> 326 master 327 </span> 328 </a> 329 </div> 330 331 <div class="select-menu-no-results">Nothing to show</div> 332 </div> 333 334 <div class="select-menu-list select-menu-tab-bucket js-select-menu-tab-bucket" data-tab-filter="tags"> 335 <div data-filterable-for="context-commitish-filter-field" data-filterable-type="substring"> 336 337 338 </div> 339 340 <div class="select-menu-no-results">Nothing to show</div> 341 </div> 342 343 </div> 344 </div> 345 </div> 346 347 <div class="btn-group right"> 348 <a href="/SasView/sasmodels/find/master" 349 class="js-show-file-finder btn btn-sm" 350 data-pjax 351 data-hotkey="t"> 352 Find file 353 </a> 354 <button aria-label="Copy file path to clipboard" class="js-zeroclipboard btn btn-sm zeroclipboard-button tooltipped tooltipped-s" data-copied-hint="Copied!" type="button">Copy path</button> 355 </div> 356 <div class="breadcrumb js-zeroclipboard-target"> 357 <span class="repo-root js-repo-root"><span itemscope="" itemtype="http://data-vocabulary.org/Breadcrumb"><a href="/SasView/sasmodels" class="" data-branch="master" data-pjax="true" itemscope="url"><span itemprop="title">sasmodels</span></a></span></span><span class="separator">/</span><span itemscope="" itemtype="http://data-vocabulary.org/Breadcrumb"><a href="/SasView/sasmodels/tree/master/sasmodels" class="" data-branch="master" data-pjax="true" itemscope="url"><span itemprop="title">sasmodels</span></a></span><span class="separator">/</span><strong class="final-path">convert.py</strong> 358 </div> 359 </div> 360 361 <include-fragment class="commit-tease" src="/SasView/sasmodels/contributors/master/sasmodels/convert.py"> 362 <div> 363 Fetching contributors… 364 </div> 365 366 <div class="commit-tease-contributors"> 367 <img alt="" class="loader-loading left" height="16" src="https://assets-cdn.github.com/images/spinners/octocat-spinner-32-EAF2F5.gif" width="16" /> 368 <span class="loader-error">Cannot retrieve contributors at this time</span> 369 </div> 370 </include-fragment> 371 <div class="file"> 372 <div class="file-header"> 373 <div class="file-actions"> 374 375 <div class="btn-group"> 376 <a href="/SasView/sasmodels/raw/master/sasmodels/convert.py" class="btn btn-sm " id="raw-url">Raw</a> 377 <a href="/SasView/sasmodels/blame/master/sasmodels/convert.py" class="btn btn-sm js-update-url-with-hash">Blame</a> 378 <a href="/SasView/sasmodels/commits/master/sasmodels/convert.py" class="btn btn-sm " rel="nofollow">History</a> 379 </div> 380 381 <a class="btn-octicon tooltipped tooltipped-nw" 382 href="https://windows.github.com" 383 aria-label="Open this file in GitHub Desktop" 384 data-ga-click="Repository, open with desktop, type:windows"> 385 <svg aria-hidden="true" class="octicon octicon-device-desktop" height="16" role="img" version="1.1" viewBox="0 0 16 16" width="16"><path d="M15 2H1c-0.55 0-1 0.45-1 1v9c0 0.55 0.45 1 1 1h5.34c-0.25 0.61-0.86 1.39-2.34 2h8c-1.48-0.61-2.09-1.39-2.34-2h5.34c0.55 0 1-0.45 1-1V3c0-0.55-0.45-1-1-1z m0 9H1V3h14v8z"></path></svg> 386 </a> 387 388 <button type="button" class="btn-octicon disabled tooltipped tooltipped-nw" 389 aria-label="You must be signed in to make or propose changes"> 390 <svg aria-hidden="true" class="octicon octicon-pencil" height="16" role="img" version="1.1" viewBox="0 0 14 16" width="14"><path d="M0 12v3h3l8-8-3-3L0 12z m3 2H1V12h1v1h1v1z m10.3-9.3l-1.3 1.3-3-3 1.3-1.3c0.39-0.39 1.02-0.39 1.41 0l1.59 1.59c0.39 0.39 0.39 1.02 0 1.41z"></path></svg> 391 </button> 392 <button type="button" class="btn-octicon btn-octicon-danger disabled tooltipped tooltipped-nw" 393 aria-label="You must be signed in to make or propose changes"> 394 <svg aria-hidden="true" class="octicon octicon-trashcan" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M10 2H8c0-0.55-0.45-1-1-1H4c-0.55 0-1 0.45-1 1H1c-0.55 0-1 0.45-1 1v1c0 0.55 0.45 1 1 1v9c0 0.55 0.45 1 1 1h7c0.55 0 1-0.45 1-1V5c0.55 0 1-0.45 1-1v-1c0-0.55-0.45-1-1-1z m-1 12H2V5h1v8h1V5h1v8h1V5h1v8h1V5h1v9z m1-10H1v-1h9v1z"></path></svg> 395 </button> 396 </div> 397 398 <div class="file-info"> 399 157 lines (141 sloc) 400 <span class="file-info-divider"></span> 401 5.04 KB 402 </div> 403 </div> 404 405 406 407 <div class="blob-wrapper data type-python"> 408 <table class="highlight tab-size js-file-line-container" data-tab-size="8"> 409 <tr> 410 <td id="L1" class="blob-num js-line-number" data-line-number="1"></td> 411 <td id="LC1" class="blob-code blob-code-inner js-file-line"><span class="pl-s"><span class="pl-pds">"""</span></span></td> 412 </tr> 413 <tr> 414 <td id="L2" class="blob-num js-line-number" data-line-number="2"></td> 415 <td id="LC2" class="blob-code blob-code-inner js-file-line"><span class="pl-s">Convert models to and from sasview.</span></td> 416 </tr> 417 <tr> 418 <td id="L3" class="blob-num js-line-number" data-line-number="3"></td> 419 <td id="LC3" class="blob-code blob-code-inner js-file-line"><span class="pl-s"><span class="pl-pds">"""</span></span></td> 420 </tr> 421 <tr> 422 <td id="L4" class="blob-num js-line-number" data-line-number="4"></td> 423 <td id="LC4" class="blob-code blob-code-inner js-file-line"><span class="pl-k">import</span> warnings</td> 424 </tr> 425 <tr> 426 <td id="L5" class="blob-num js-line-number" data-line-number="5"></td> 427 <td id="LC5" class="blob-code blob-code-inner js-file-line"> 428 </td> 429 </tr> 430 <tr> 431 <td id="L6" class="blob-num js-line-number" data-line-number="6"></td> 432 <td id="LC6" class="blob-code blob-code-inner js-file-line"><span class="pl-c"># List of models which SasView versions don't contain the explicit 'scale' argument.</span></td> 433 </tr> 434 <tr> 435 <td id="L7" class="blob-num js-line-number" data-line-number="7"></td> 436 <td id="LC7" class="blob-code blob-code-inner js-file-line"><span class="pl-c"># When converting such a model, please update this list.</span></td> 437 </tr> 438 <tr> 439 <td id="L8" class="blob-num js-line-number" data-line-number="8"></td> 440 <td id="LC8" class="blob-code blob-code-inner js-file-line"><span class="pl-c1">MODELS_WITHOUT_SCALE</span> <span class="pl-k">=</span> [</td> 441 </tr> 442 <tr> 443 <td id="L9" class="blob-num js-line-number" data-line-number="9"></td> 444 <td id="LC9" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>teubner_strey<span class="pl-pds">'</span></span>,</td> 445 </tr> 446 <tr> 447 <td id="L10" class="blob-num js-line-number" data-line-number="10"></td> 448 <td id="LC10" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>broad_peak<span class="pl-pds">'</span></span>,</td> 449 </tr> 450 <tr> 451 <td id="L11" class="blob-num js-line-number" data-line-number="11"></td> 452 <td id="LC11" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>two_lorentzian<span class="pl-pds">'</span></span>,</td> 453 </tr> 454 <tr> 455 <td id="L12" class="blob-num js-line-number" data-line-number="12"></td> 456 <td id="LC12" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>gel_fit<span class="pl-pds">'</span></span>,</td> 457 </tr> 458 <tr> 459 <td id="L13" class="blob-num js-line-number" data-line-number="13"></td> 460 <td id="LC13" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>gauss_lorentz_gel<span class="pl-pds">'</span></span>,</td> 461 </tr> 462 <tr> 463 <td id="L14" class="blob-num js-line-number" data-line-number="14"></td> 464 <td id="LC14" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>be_polyelectrolyte<span class="pl-pds">'</span></span>,</td> 465 </tr> 466 <tr> 467 <td id="L15" class="blob-num js-line-number" data-line-number="15"></td> 468 <td id="LC15" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>correlation_length<span class="pl-pds">'</span></span>,</td> 469 </tr> 470 <tr> 471 <td id="L16" class="blob-num js-line-number" data-line-number="16"></td> 472 <td id="LC16" class="blob-code blob-code-inner js-file-line">]</td> 473 </tr> 474 <tr> 475 <td id="L17" class="blob-num js-line-number" data-line-number="17"></td> 476 <td id="LC17" class="blob-code blob-code-inner js-file-line"> 477 </td> 478 </tr> 479 <tr> 480 <td id="L18" class="blob-num js-line-number" data-line-number="18"></td> 481 <td id="LC18" class="blob-code blob-code-inner js-file-line"><span class="pl-c"># List of models which SasView versions don't contain the explicit 'background' argument.</span></td> 482 </tr> 483 <tr> 484 <td id="L19" class="blob-num js-line-number" data-line-number="19"></td> 485 <td id="LC19" class="blob-code blob-code-inner js-file-line"><span class="pl-c"># When converting such a model, please update this list.</span></td> 486 </tr> 487 <tr> 488 <td id="L20" class="blob-num js-line-number" data-line-number="20"></td> 489 <td id="LC20" class="blob-code blob-code-inner js-file-line"><span class="pl-c1">MODELS_WITHOUT_BACKGROUND</span> <span class="pl-k">=</span> [</td> 490 </tr> 491 <tr> 492 <td id="L21" class="blob-num js-line-number" data-line-number="21"></td> 493 <td id="LC21" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">'</span>guinier<span class="pl-pds">'</span></span>,</td> 494 </tr> 495 <tr> 496 <td id="L22" class="blob-num js-line-number" data-line-number="22"></td> 497 <td id="LC22" class="blob-code blob-code-inner js-file-line">]</td> 498 </tr> 499 <tr> 500 <td id="L23" class="blob-num js-line-number" data-line-number="23"></td> 501 <td id="LC23" class="blob-code blob-code-inner js-file-line"> 502 </td> 503 </tr> 504 <tr> 505 <td id="L24" class="blob-num js-line-number" data-line-number="24"></td> 506 <td id="LC24" class="blob-code blob-code-inner js-file-line"><span class="pl-c1">PD_DOT</span> <span class="pl-k">=</span> [</td> 507 </tr> 508 <tr> 509 <td id="L25" class="blob-num js-line-number" data-line-number="25"></td> 510 <td id="LC25" class="blob-code blob-code-inner js-file-line"> (<span class="pl-s"><span class="pl-pds">"</span><span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span><span class="pl-pds">"</span></span>),</td> 511 </tr> 512 <tr> 513 <td id="L26" class="blob-num js-line-number" data-line-number="26"></td> 514 <td id="LC26" class="blob-code blob-code-inner js-file-line"> (<span class="pl-s"><span class="pl-pds">"</span>_pd<span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span>.width<span class="pl-pds">"</span></span>),</td> 515 </tr> 516 <tr> 517 <td id="L27" class="blob-num js-line-number" data-line-number="27"></td> 518 <td id="LC27" class="blob-code blob-code-inner js-file-line"> (<span class="pl-s"><span class="pl-pds">"</span>_pd_n<span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span>.npts<span class="pl-pds">"</span></span>),</td> 519 </tr> 520 <tr> 521 <td id="L28" class="blob-num js-line-number" data-line-number="28"></td> 522 <td id="LC28" class="blob-code blob-code-inner js-file-line"> (<span class="pl-s"><span class="pl-pds">"</span>_pd_nsigma<span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span>.nsigmas<span class="pl-pds">"</span></span>),</td> 523 </tr> 524 <tr> 525 <td id="L29" class="blob-num js-line-number" data-line-number="29"></td> 526 <td id="LC29" class="blob-code blob-code-inner js-file-line"> (<span class="pl-s"><span class="pl-pds">"</span>_pd_type<span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span>.type<span class="pl-pds">"</span></span>),</td> 527 </tr> 528 <tr> 529 <td id="L30" class="blob-num js-line-number" data-line-number="30"></td> 530 <td id="LC30" class="blob-code blob-code-inner js-file-line"> ]</td> 531 </tr> 532 <tr> 533 <td id="L31" class="blob-num js-line-number" data-line-number="31"></td> 534 <td id="LC31" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">_convert_pars</span>(<span class="pl-smi">pars</span>, <span class="pl-smi">mapping</span>):</td> 535 </tr> 536 <tr> 537 <td id="L32" class="blob-num js-line-number" data-line-number="32"></td> 538 <td id="LC32" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 539 </tr> 540 <tr> 541 <td id="L33" class="blob-num js-line-number" data-line-number="33"></td> 542 <td id="LC33" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> Rename the parameters and any associated polydispersity attributes.</span></td> 543 </tr> 544 <tr> 545 <td id="L34" class="blob-num js-line-number" data-line-number="34"></td> 546 <td id="LC34" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 547 </tr> 548 <tr> 549 <td id="L35" class="blob-num js-line-number" data-line-number="35"></td> 550 <td id="LC35" class="blob-code blob-code-inner js-file-line"> newpars <span class="pl-k">=</span> pars.copy()</td> 551 </tr> 552 <tr> 553 <td id="L36" class="blob-num js-line-number" data-line-number="36"></td> 554 <td id="LC36" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> new, old <span class="pl-k">in</span> mapping.items():</td> 555 </tr> 556 <tr> 557 <td id="L37" class="blob-num js-line-number" data-line-number="37"></td> 558 <td id="LC37" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> old <span class="pl-k">==</span> new: <span class="pl-k">continue</span></td> 559 </tr> 560 <tr> 561 <td id="L38" class="blob-num js-line-number" data-line-number="38"></td> 562 <td id="LC38" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> pd, dot <span class="pl-k">in</span> <span class="pl-c1">PD_DOT</span>:</td> 563 </tr> 564 <tr> 565 <td id="L39" class="blob-num js-line-number" data-line-number="39"></td> 566 <td id="LC39" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> old<span class="pl-k">+</span>dot <span class="pl-k">in</span> newpars:</td> 567 </tr> 568 <tr> 569 <td id="L40" class="blob-num js-line-number" data-line-number="40"></td> 570 <td id="LC40" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> new <span class="pl-k">is</span> <span class="pl-k">not</span> <span class="pl-c1">None</span>:</td> 571 </tr> 572 <tr> 573 <td id="L41" class="blob-num js-line-number" data-line-number="41"></td> 574 <td id="LC41" class="blob-code blob-code-inner js-file-line"> newpars[new<span class="pl-k">+</span>pd] <span class="pl-k">=</span> pars[old<span class="pl-k">+</span>dot]</td> 575 </tr> 576 <tr> 577 <td id="L42" class="blob-num js-line-number" data-line-number="42"></td> 578 <td id="LC42" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">del</span> newpars[old<span class="pl-k">+</span>dot]</td> 579 </tr> 580 <tr> 581 <td id="L43" class="blob-num js-line-number" data-line-number="43"></td> 582 <td id="LC43" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">return</span> newpars</td> 583 </tr> 584 <tr> 585 <td id="L44" class="blob-num js-line-number" data-line-number="44"></td> 586 <td id="LC44" class="blob-code blob-code-inner js-file-line"> 587 </td> 588 </tr> 589 <tr> 590 <td id="L45" class="blob-num js-line-number" data-line-number="45"></td> 591 <td id="LC45" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">_rescale_sld</span>(<span class="pl-smi">pars</span>):</td> 592 </tr> 593 <tr> 594 <td id="L46" class="blob-num js-line-number" data-line-number="46"></td> 595 <td id="LC46" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 596 </tr> 597 <tr> 598 <td id="L47" class="blob-num js-line-number" data-line-number="47"></td> 599 <td id="LC47" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> rescale all sld parameters in the new model definition by 1e6 so the</span></td> 600 </tr> 601 <tr> 602 <td id="L48" class="blob-num js-line-number" data-line-number="48"></td> 603 <td id="LC48" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> numbers are nicer. Relies on the fact that all sld parameters in the</span></td> 604 </tr> 605 <tr> 606 <td id="L49" class="blob-num js-line-number" data-line-number="49"></td> 607 <td id="LC49" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> new model definition end with sld.</span></td> 608 </tr> 609 <tr> 610 <td id="L50" class="blob-num js-line-number" data-line-number="50"></td> 611 <td id="LC50" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 612 </tr> 613 <tr> 614 <td id="L51" class="blob-num js-line-number" data-line-number="51"></td> 615 <td id="LC51" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">return</span> <span class="pl-c1">dict</span>((p, (v<span class="pl-k">*</span><span class="pl-c1">1e6</span> <span class="pl-k">if</span> p.endswith(<span class="pl-s"><span class="pl-pds">'</span>sld<span class="pl-pds">'</span></span>) <span class="pl-k">else</span> v))</td> 616 </tr> 617 <tr> 618 <td id="L52" class="blob-num js-line-number" data-line-number="52"></td> 619 <td id="LC52" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> p, v <span class="pl-k">in</span> pars.items())</td> 620 </tr> 621 <tr> 622 <td id="L53" class="blob-num js-line-number" data-line-number="53"></td> 623 <td id="LC53" class="blob-code blob-code-inner js-file-line"> 624 </td> 625 </tr> 626 <tr> 627 <td id="L54" class="blob-num js-line-number" data-line-number="54"></td> 628 <td id="LC54" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">convert_model</span>(<span class="pl-smi">name</span>, <span class="pl-smi">pars</span>):</td> 629 </tr> 630 <tr> 631 <td id="L55" class="blob-num js-line-number" data-line-number="55"></td> 632 <td id="LC55" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 633 </tr> 634 <tr> 635 <td id="L56" class="blob-num js-line-number" data-line-number="56"></td> 636 <td id="LC56" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> Convert model from old style parameter names to new style.</span></td> 637 </tr> 638 <tr> 639 <td id="L57" class="blob-num js-line-number" data-line-number="57"></td> 640 <td id="LC57" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 641 </tr> 642 <tr> 643 <td id="L58" class="blob-num js-line-number" data-line-number="58"></td> 644 <td id="LC58" class="blob-code blob-code-inner js-file-line"> _, _ <span class="pl-k">=</span> name, pars <span class="pl-c"># lint</span></td> 645 </tr> 646 <tr> 647 <td id="L59" class="blob-num js-line-number" data-line-number="59"></td> 648 <td id="LC59" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">raise</span> <span class="pl-c1">NotImplementedError</span></td> 649 </tr> 650 <tr> 651 <td id="L60" class="blob-num js-line-number" data-line-number="60"></td> 652 <td id="LC60" class="blob-code blob-code-inner js-file-line"> <span class="pl-c"># need to load all new models in order to determine old=>new</span></td> 653 </tr> 654 <tr> 655 <td id="L61" class="blob-num js-line-number" data-line-number="61"></td> 656 <td id="LC61" class="blob-code blob-code-inner js-file-line"> <span class="pl-c"># model name mapping</span></td> 657 </tr> 658 <tr> 659 <td id="L62" class="blob-num js-line-number" data-line-number="62"></td> 660 <td id="LC62" class="blob-code blob-code-inner js-file-line"> 661 </td> 662 </tr> 663 <tr> 664 <td id="L63" class="blob-num js-line-number" data-line-number="63"></td> 665 <td id="LC63" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">_unscale_sld</span>(<span class="pl-smi">pars</span>):</td> 666 </tr> 667 <tr> 668 <td id="L64" class="blob-num js-line-number" data-line-number="64"></td> 669 <td id="LC64" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 670 </tr> 671 <tr> 672 <td id="L65" class="blob-num js-line-number" data-line-number="65"></td> 673 <td id="LC65" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> rescale all sld parameters in the new model definition by 1e6 so the</span></td> 674 </tr> 675 <tr> 676 <td id="L66" class="blob-num js-line-number" data-line-number="66"></td> 677 <td id="LC66" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> numbers are nicer. Relies on the fact that all sld parameters in the</span></td> 678 </tr> 679 <tr> 680 <td id="L67" class="blob-num js-line-number" data-line-number="67"></td> 681 <td id="LC67" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> new model definition end with sld.</span></td> 682 </tr> 683 <tr> 684 <td id="L68" class="blob-num js-line-number" data-line-number="68"></td> 685 <td id="LC68" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 686 </tr> 687 <tr> 688 <td id="L69" class="blob-num js-line-number" data-line-number="69"></td> 689 <td id="LC69" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">return</span> <span class="pl-c1">dict</span>((p, (v<span class="pl-k">*</span><span class="pl-c1">1e-6</span> <span class="pl-k">if</span> p.endswith(<span class="pl-s"><span class="pl-pds">'</span>sld<span class="pl-pds">'</span></span>) <span class="pl-k">else</span> v))</td> 690 </tr> 691 <tr> 692 <td id="L70" class="blob-num js-line-number" data-line-number="70"></td> 693 <td id="LC70" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> p, v <span class="pl-k">in</span> pars.items())</td> 694 </tr> 695 <tr> 696 <td id="L71" class="blob-num js-line-number" data-line-number="71"></td> 697 <td id="LC71" class="blob-code blob-code-inner js-file-line"> 698 </td> 699 </tr> 700 <tr> 701 <td id="L72" class="blob-num js-line-number" data-line-number="72"></td> 702 <td id="LC72" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">_remove_pd</span>(<span class="pl-smi">pars</span>, <span class="pl-smi">key</span>, <span class="pl-smi">name</span>):</td> 703 </tr> 704 <tr> 705 <td id="L73" class="blob-num js-line-number" data-line-number="73"></td> 706 <td id="LC73" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 707 </tr> 708 <tr> 709 <td id="L74" class="blob-num js-line-number" data-line-number="74"></td> 710 <td id="LC74" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> Remove polydispersity from the parameter list.</span></td> 711 </tr> 712 <tr> 713 <td id="L75" class="blob-num js-line-number" data-line-number="75"></td> 714 <td id="LC75" class="blob-code blob-code-inner js-file-line"><span class="pl-s"></span></td> 715 </tr> 716 <tr> 717 <td id="L76" class="blob-num js-line-number" data-line-number="76"></td> 718 <td id="LC76" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> Note: operates in place</span></td> 719 </tr> 720 <tr> 721 <td id="L77" class="blob-num js-line-number" data-line-number="77"></td> 722 <td id="LC77" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 723 </tr> 724 <tr> 725 <td id="L78" class="blob-num js-line-number" data-line-number="78"></td> 726 <td id="LC78" class="blob-code blob-code-inner js-file-line"> <span class="pl-c"># Bumps style parameter names</span></td> 727 </tr> 728 <tr> 729 <td id="L79" class="blob-num js-line-number" data-line-number="79"></td> 730 <td id="LC79" class="blob-code blob-code-inner js-file-line"> pd <span class="pl-k">=</span> pars.pop(key<span class="pl-k">+</span><span class="pl-s"><span class="pl-pds">"</span>.width<span class="pl-pds">"</span></span>, <span class="pl-c1">0.0</span>)</td> 731 </tr> 732 <tr> 733 <td id="L80" class="blob-num js-line-number" data-line-number="80"></td> 734 <td id="LC80" class="blob-code blob-code-inner js-file-line"> pd_n <span class="pl-k">=</span> pars.pop(key<span class="pl-k">+</span><span class="pl-s"><span class="pl-pds">"</span>.npts<span class="pl-pds">"</span></span>, <span class="pl-c1">0</span>)</td> 735 </tr> 736 <tr> 737 <td id="L81" class="blob-num js-line-number" data-line-number="81"></td> 738 <td id="LC81" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> pd <span class="pl-k">!=</span> <span class="pl-c1">0.0</span> <span class="pl-k">and</span> pd_n <span class="pl-k">!=</span> <span class="pl-c1">0</span>:</td> 739 </tr> 740 <tr> 741 <td id="L82" class="blob-num js-line-number" data-line-number="82"></td> 742 <td id="LC82" class="blob-code blob-code-inner js-file-line"> warnings.warn(<span class="pl-s"><span class="pl-pds">"</span>parameter <span class="pl-c1">%s</span> not polydisperse in sasview <span class="pl-c1">%s</span><span class="pl-pds">"</span></span><span class="pl-k">%</span>(key, name))</td> 743 </tr> 744 <tr> 745 <td id="L83" class="blob-num js-line-number" data-line-number="83"></td> 746 <td id="LC83" class="blob-code blob-code-inner js-file-line"> pars.pop(key<span class="pl-k">+</span><span class="pl-s"><span class="pl-pds">"</span>.nsigmas<span class="pl-pds">"</span></span>, <span class="pl-c1">None</span>)</td> 747 </tr> 748 <tr> 749 <td id="L84" class="blob-num js-line-number" data-line-number="84"></td> 750 <td id="LC84" class="blob-code blob-code-inner js-file-line"> pars.pop(key<span class="pl-k">+</span><span class="pl-s"><span class="pl-pds">"</span>.type<span class="pl-pds">"</span></span>, <span class="pl-c1">None</span>)</td> 751 </tr> 752 <tr> 753 <td id="L85" class="blob-num js-line-number" data-line-number="85"></td> 754 <td id="LC85" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">return</span> pars</td> 755 </tr> 756 <tr> 757 <td id="L86" class="blob-num js-line-number" data-line-number="86"></td> 758 <td id="LC86" class="blob-code blob-code-inner js-file-line"> 759 </td> 760 </tr> 761 <tr> 762 <td id="L87" class="blob-num js-line-number" data-line-number="87"></td> 763 <td id="LC87" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">_revert_pars</span>(<span class="pl-smi">pars</span>, <span class="pl-smi">mapping</span>):</td> 764 </tr> 765 <tr> 766 <td id="L88" class="blob-num js-line-number" data-line-number="88"></td> 767 <td id="LC88" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 768 </tr> 769 <tr> 770 <td id="L89" class="blob-num js-line-number" data-line-number="89"></td> 771 <td id="LC89" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> Rename the parameters and any associated polydispersity attributes.</span></td> 772 </tr> 773 <tr> 774 <td id="L90" class="blob-num js-line-number" data-line-number="90"></td> 775 <td id="LC90" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 776 </tr> 777 <tr> 778 <td id="L91" class="blob-num js-line-number" data-line-number="91"></td> 779 <td id="LC91" class="blob-code blob-code-inner js-file-line"> newpars <span class="pl-k">=</span> pars.copy()</td> 780 </tr> 781 <tr> 782 <td id="L92" class="blob-num js-line-number" data-line-number="92"></td> 783 <td id="LC92" class="blob-code blob-code-inner js-file-line"> 784 </td> 785 </tr> 786 <tr> 787 <td id="L93" class="blob-num js-line-number" data-line-number="93"></td> 788 <td id="LC93" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> new, old <span class="pl-k">in</span> mapping.items():</td> 789 </tr> 790 <tr> 791 <td id="L94" class="blob-num js-line-number" data-line-number="94"></td> 792 <td id="LC94" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> pd, dot <span class="pl-k">in</span> <span class="pl-c1">PD_DOT</span>:</td> 793 </tr> 794 <tr> 795 <td id="L95" class="blob-num js-line-number" data-line-number="95"></td> 796 <td id="LC95" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> old <span class="pl-k">and</span> old<span class="pl-k">+</span>pd <span class="pl-k">==</span> new<span class="pl-k">+</span>dot:</td> 797 </tr> 798 <tr> 799 <td id="L96" class="blob-num js-line-number" data-line-number="96"></td> 800 <td id="LC96" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">continue</span></td> 801 </tr> 802 <tr> 803 <td id="L97" class="blob-num js-line-number" data-line-number="97"></td> 804 <td id="LC97" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> new<span class="pl-k">+</span>pd <span class="pl-k">in</span> newpars:</td> 805 </tr> 806 <tr> 807 <td id="L98" class="blob-num js-line-number" data-line-number="98"></td> 808 <td id="LC98" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> old <span class="pl-k">is</span> <span class="pl-k">not</span> <span class="pl-c1">None</span>:</td> 809 </tr> 810 <tr> 811 <td id="L99" class="blob-num js-line-number" data-line-number="99"></td> 812 <td id="LC99" class="blob-code blob-code-inner js-file-line"> newpars[old<span class="pl-k">+</span>dot] <span class="pl-k">=</span> pars[new<span class="pl-k">+</span>pd]</td> 813 </tr> 814 <tr> 815 <td id="L100" class="blob-num js-line-number" data-line-number="100"></td> 816 <td id="LC100" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">del</span> newpars[new<span class="pl-k">+</span>pd]</td> 817 </tr> 818 <tr> 819 <td id="L101" class="blob-num js-line-number" data-line-number="101"></td> 820 <td id="LC101" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> k <span class="pl-k">in</span> <span class="pl-c1">list</span>(newpars.keys()):</td> 821 </tr> 822 <tr> 823 <td id="L102" class="blob-num js-line-number" data-line-number="102"></td> 824 <td id="LC102" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> pd, dot <span class="pl-k">in</span> <span class="pl-c1">PD_DOT</span>[<span class="pl-c1">1</span>:]: <span class="pl-c"># skip "" => ""</span></td> 825 </tr> 826 <tr> 827 <td id="L103" class="blob-num js-line-number" data-line-number="103"></td> 828 <td id="LC103" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> k.endswith(pd):</td> 829 </tr> 830 <tr> 831 <td id="L104" class="blob-num js-line-number" data-line-number="104"></td> 832 <td id="LC104" class="blob-code blob-code-inner js-file-line"> newpars[k[:<span class="pl-k">-</span><span class="pl-c1">len</span>(pd)]<span class="pl-k">+</span>dot] <span class="pl-k">=</span> newpars[k]</td> 833 </tr> 834 <tr> 835 <td id="L105" class="blob-num js-line-number" data-line-number="105"></td> 836 <td id="LC105" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">del</span> newpars[k]</td> 837 </tr> 838 <tr> 839 <td id="L106" class="blob-num js-line-number" data-line-number="106"></td> 840 <td id="LC106" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">return</span> newpars</td> 841 </tr> 842 <tr> 843 <td id="L107" class="blob-num js-line-number" data-line-number="107"></td> 844 <td id="LC107" class="blob-code blob-code-inner js-file-line"> 845 </td> 846 </tr> 847 <tr> 848 <td id="L108" class="blob-num js-line-number" data-line-number="108"></td> 849 <td id="LC108" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">revert_model</span>(<span class="pl-smi">model_definition</span>, <span class="pl-smi">pars</span>):</td> 850 </tr> 851 <tr> 852 <td id="L109" class="blob-num js-line-number" data-line-number="109"></td> 853 <td id="LC109" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 854 </tr> 855 <tr> 856 <td id="L110" class="blob-num js-line-number" data-line-number="110"></td> 857 <td id="LC110" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> Convert model from new style parameter names to old style.</span></td> 858 </tr> 859 <tr> 860 <td id="L111" class="blob-num js-line-number" data-line-number="111"></td> 861 <td id="LC111" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 862 </tr> 863 <tr> 864 <td id="L112" class="blob-num js-line-number" data-line-number="112"></td> 865 <td id="LC112" class="blob-code blob-code-inner js-file-line"> mapping <span class="pl-k">=</span> model_definition.oldpars</td> 866 </tr> 867 <tr> 868 <td id="L113" class="blob-num js-line-number" data-line-number="113"></td> 869 <td id="LC113" class="blob-code blob-code-inner js-file-line"> oldname <span class="pl-k">=</span> model_definition.oldname</td> 870 </tr> 871 <tr> 872 <td id="L114" class="blob-num js-line-number" data-line-number="114"></td> 873 <td id="LC114" class="blob-code blob-code-inner js-file-line"> oldpars <span class="pl-k">=</span> _revert_pars(_unscale_sld(pars), mapping)</td> 874 </tr> 875 <tr> 876 <td id="L115" class="blob-num js-line-number" data-line-number="115"></td> 877 <td id="LC115" class="blob-code blob-code-inner js-file-line"> 878 </td> 879 </tr> 880 <tr> 881 <td id="L116" class="blob-num js-line-number" data-line-number="116"></td> 882 <td id="LC116" class="blob-code blob-code-inner js-file-line"> <span class="pl-c"># Note: update compare.constrain_pars to match</span></td> 883 </tr> 884 <tr> 885 <td id="L117" class="blob-num js-line-number" data-line-number="117"></td> 886 <td id="LC117" class="blob-code blob-code-inner js-file-line"> name <span class="pl-k">=</span> model_definition.name</td> 887 </tr> 888 <tr> 889 <td id="L118" class="blob-num js-line-number" data-line-number="118"></td> 890 <td id="LC118" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> name <span class="pl-k">in</span> <span class="pl-c1">MODELS_WITHOUT_SCALE</span>:</td> 891 </tr> 892 <tr> 893 <td id="L119" class="blob-num js-line-number" data-line-number="119"></td> 894 <td id="LC119" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> oldpars.pop(<span class="pl-s"><span class="pl-pds">'</span>scale<span class="pl-pds">'</span></span>, <span class="pl-c1">1.0</span>) <span class="pl-k">!=</span> <span class="pl-c1">1.0</span>:</td> 895 </tr> 896 <tr> 897 <td id="L120" class="blob-num js-line-number" data-line-number="120"></td> 898 <td id="LC120" class="blob-code blob-code-inner js-file-line"> warnings.warn(<span class="pl-s"><span class="pl-pds">"</span>parameter scale not used in sasview <span class="pl-c1">%s</span><span class="pl-pds">"</span></span><span class="pl-k">%</span>name)</td> 899 </tr> 900 <tr> 901 <td id="L121" class="blob-num js-line-number" data-line-number="121"></td> 902 <td id="LC121" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> name <span class="pl-k">in</span> <span class="pl-c1">MODELS_WITHOUT_BACKGROUND</span>:</td> 903 </tr> 904 <tr> 905 <td id="L122" class="blob-num js-line-number" data-line-number="122"></td> 906 <td id="LC122" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> oldpars.pop(<span class="pl-s"><span class="pl-pds">'</span>background<span class="pl-pds">'</span></span>, <span class="pl-c1">0.0</span>) <span class="pl-k">!=</span> <span class="pl-c1">0.0</span>:</td> 907 </tr> 908 <tr> 909 <td id="L123" class="blob-num js-line-number" data-line-number="123"></td> 910 <td id="LC123" class="blob-code blob-code-inner js-file-line"> warnings.warn(<span class="pl-s"><span class="pl-pds">"</span>parameter background not used in sasview <span class="pl-c1">%s</span><span class="pl-pds">"</span></span><span class="pl-k">%</span>name)</td> 911 </tr> 912 <tr> 913 <td id="L124" class="blob-num js-line-number" data-line-number="124"></td> 914 <td id="LC124" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> <span class="pl-c1">getattr</span>(model_definition, <span class="pl-s"><span class="pl-pds">'</span>category<span class="pl-pds">'</span></span>, <span class="pl-c1">None</span>) <span class="pl-k">==</span> <span class="pl-s"><span class="pl-pds">'</span>structure-factor<span class="pl-pds">'</span></span>:</td> 915 </tr> 916 <tr> 917 <td id="L125" class="blob-num js-line-number" data-line-number="125"></td> 918 <td id="LC125" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> oldpars.pop(<span class="pl-s"><span class="pl-pds">'</span>scale<span class="pl-pds">'</span></span>, <span class="pl-c1">1.0</span>) <span class="pl-k">!=</span> <span class="pl-c1">1.0</span>:</td> 919 </tr> 920 <tr> 921 <td id="L126" class="blob-num js-line-number" data-line-number="126"></td> 922 <td id="LC126" class="blob-code blob-code-inner js-file-line"> warnings.warn(<span class="pl-s"><span class="pl-pds">"</span>parameter scale not used in sasview <span class="pl-c1">%s</span><span class="pl-pds">"</span></span><span class="pl-k">%</span>name)</td> 923 </tr> 924 <tr> 925 <td id="L127" class="blob-num js-line-number" data-line-number="127"></td> 926 <td id="LC127" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> oldpars.pop(<span class="pl-s"><span class="pl-pds">'</span>background<span class="pl-pds">'</span></span>, <span class="pl-c1">0.0</span>) <span class="pl-k">!=</span> <span class="pl-c1">0.0</span>:</td> 927 </tr> 928 <tr> 929 <td id="L128" class="blob-num js-line-number" data-line-number="128"></td> 930 <td id="LC128" class="blob-code blob-code-inner js-file-line"> warnings.warn(<span class="pl-s"><span class="pl-pds">"</span>parameter background not used in sasview <span class="pl-c1">%s</span><span class="pl-pds">"</span></span><span class="pl-k">%</span>name)</td> 931 </tr> 932 <tr> 933 <td id="L129" class="blob-num js-line-number" data-line-number="129"></td> 934 <td id="LC129" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> name <span class="pl-k">==</span> <span class="pl-s"><span class="pl-pds">'</span>pearl_necklace<span class="pl-pds">'</span></span>:</td> 935 </tr> 936 <tr> 937 <td id="L130" class="blob-num js-line-number" data-line-number="130"></td> 938 <td id="LC130" class="blob-code blob-code-inner js-file-line"> _remove_pd(oldpars, <span class="pl-s"><span class="pl-pds">'</span>num_pearls<span class="pl-pds">'</span></span>, name)</td> 939 </tr> 940 <tr> 941 <td id="L131" class="blob-num js-line-number" data-line-number="131"></td> 942 <td id="LC131" class="blob-code blob-code-inner js-file-line"> _remove_pd(oldpars, <span class="pl-s"><span class="pl-pds">'</span>thick_string<span class="pl-pds">'</span></span>, name)</td> 943 </tr> 944 <tr> 945 <td id="L132" class="blob-num js-line-number" data-line-number="132"></td> 946 <td id="LC132" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> name <span class="pl-k">==</span> <span class="pl-s"><span class="pl-pds">'</span>rpa<span class="pl-pds">'</span></span>:</td> 947 </tr> 948 <tr> 949 <td id="L133" class="blob-num js-line-number" data-line-number="133"></td> 950 <td id="LC133" class="blob-code blob-code-inner js-file-line"> <span class="pl-c"># convert scattering lengths from femtometers to centimeters</span></td> 951 </tr> 952 <tr> 953 <td id="L134" class="blob-num js-line-number" data-line-number="134"></td> 954 <td id="LC134" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">for</span> p <span class="pl-k">in</span> <span class="pl-s"><span class="pl-pds">"</span>La<span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span>Lb<span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span>Lc<span class="pl-pds">"</span></span>, <span class="pl-s"><span class="pl-pds">"</span>Ld<span class="pl-pds">"</span></span>:</td> 955 </tr> 956 <tr> 957 <td id="L135" class="blob-num js-line-number" data-line-number="135"></td> 958 <td id="LC135" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> p <span class="pl-k">in</span> oldpars: oldpars[p] <span class="pl-k">*=</span> <span class="pl-c1">1e-13</span></td> 959 </tr> 960 <tr> 961 <td id="L136" class="blob-num js-line-number" data-line-number="136"></td> 962 <td id="LC136" class="blob-code blob-code-inner js-file-line"> 963 </td> 964 </tr> 965 <tr> 966 <td id="L137" class="blob-num js-line-number" data-line-number="137"></td> 967 <td id="LC137" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">return</span> oldname, oldpars</td> 968 </tr> 969 <tr> 970 <td id="L138" class="blob-num js-line-number" data-line-number="138"></td> 971 <td id="LC138" class="blob-code blob-code-inner js-file-line"> 972 </td> 973 </tr> 974 <tr> 975 <td id="L139" class="blob-num js-line-number" data-line-number="139"></td> 976 <td id="LC139" class="blob-code blob-code-inner js-file-line"><span class="pl-k">def</span> <span class="pl-en">constrain_new_to_old</span>(<span class="pl-smi">model_definition</span>, <span class="pl-smi">pars</span>):</td> 977 </tr> 978 <tr> 979 <td id="L140" class="blob-num js-line-number" data-line-number="140"></td> 980 <td id="LC140" class="blob-code blob-code-inner js-file-line"> <span class="pl-s"><span class="pl-pds">"""</span></span></td> 981 </tr> 982 <tr> 983 <td id="L141" class="blob-num js-line-number" data-line-number="141"></td> 984 <td id="LC141" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> Restrict parameter values to those that will match sasview.</span></td> 985 </tr> 986 <tr> 987 <td id="L142" class="blob-num js-line-number" data-line-number="142"></td> 988 <td id="LC142" class="blob-code blob-code-inner js-file-line"><span class="pl-s"> <span class="pl-pds">"""</span></span></td> 989 </tr> 990 <tr> 991 <td id="L143" class="blob-num js-line-number" data-line-number="143"></td> 992 <td id="LC143" class="blob-code blob-code-inner js-file-line"> <span class="pl-c"># Note: update convert.revert_model to match</span></td> 993 </tr> 994 <tr> 995 <td id="L144" class="blob-num js-line-number" data-line-number="144"></td> 996 <td id="LC144" class="blob-code blob-code-inner js-file-line"> name <span class="pl-k">=</span> model_definition.name</td> 997 </tr> 998 <tr> 999 <td id="L145" class="blob-num js-line-number" data-line-number="145"></td> 1000 <td id="LC145" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">if</span> name <span class="pl-k">in</span> <span class="pl-c1">MODELS_WITHOUT_SCALE</span>:</td> 1001 </tr> 1002 <tr> 1003 <td id="L146" class="blob-num js-line-number" data-line-number="146"></td> 1004 <td id="LC146" class="blob-code blob-code-inner js-file-line"> pars[<span class="pl-s"><span class="pl-pds">'</span>scale<span class="pl-pds">'</span></span>] <span class="pl-k">=</span> <span class="pl-c1">1</span></td> 1005 </tr> 1006 <tr> 1007 <td id="L147" class="blob-num js-line-number" data-line-number="147"></td> 1008 <td id="LC147" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> name <span class="pl-k">in</span> <span class="pl-c1">MODELS_WITHOUT_BACKGROUND</span>:</td> 1009 </tr> 1010 <tr> 1011 <td id="L148" class="blob-num js-line-number" data-line-number="148"></td> 1012 <td id="LC148" class="blob-code blob-code-inner js-file-line"> pars[<span class="pl-s"><span class="pl-pds">'</span>background<span class="pl-pds">'</span></span>] <span class="pl-k">=</span> <span class="pl-c1">0</span></td> 1013 </tr> 1014 <tr> 1015 <td id="L149" class="blob-num js-line-number" data-line-number="149"></td> 1016 <td id="LC149" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> name <span class="pl-k">==</span> <span class="pl-s"><span class="pl-pds">'</span>pearl_necklace<span class="pl-pds">'</span></span>:</td> 1017 </tr> 1018 <tr> 1019 <td id="L150" class="blob-num js-line-number" data-line-number="150"></td> 1020 <td id="LC150" class="blob-code blob-code-inner js-file-line"> pars[<span class="pl-s"><span class="pl-pds">'</span>string_thickness_pd_n<span class="pl-pds">'</span></span>] <span class="pl-k">=</span> <span class="pl-c1">0</span></td> 1021 </tr> 1022 <tr> 1023 <td id="L151" class="blob-num js-line-number" data-line-number="151"></td> 1024 <td id="LC151" class="blob-code blob-code-inner js-file-line"> pars[<span class="pl-s"><span class="pl-pds">'</span>number_of_pearls_pd_n<span class="pl-pds">'</span></span>] <span class="pl-k">=</span> <span class="pl-c1">0</span></td> 1025 </tr> 1026 <tr> 1027 <td id="L152" class="blob-num js-line-number" data-line-number="152"></td> 1028 <td id="LC152" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> name <span class="pl-k">==</span> <span class="pl-s"><span class="pl-pds">'</span>rpa<span class="pl-pds">'</span></span>:</td> 1029 </tr> 1030 <tr> 1031 <td id="L153" class="blob-num js-line-number" data-line-number="153"></td> 1032 <td id="LC153" class="blob-code blob-code-inner js-file-line"> pars[<span class="pl-s"><span class="pl-pds">'</span>case_num<span class="pl-pds">'</span></span>] <span class="pl-k">=</span> <span class="pl-c1">int</span>(pars[<span class="pl-s"><span class="pl-pds">'</span>case_num<span class="pl-pds">'</span></span>])</td> 1033 </tr> 1034 <tr> 1035 <td id="L154" class="blob-num js-line-number" data-line-number="154"></td> 1036 <td id="LC154" class="blob-code blob-code-inner js-file-line"> <span class="pl-k">elif</span> <span class="pl-c1">getattr</span>(model_definition, <span class="pl-s"><span class="pl-pds">'</span>category<span class="pl-pds">'</span></span>, <span class="pl-c1">None</span>) <span class="pl-k">==</span> <span class="pl-s"><span class="pl-pds">'</span>structure-factor<span class="pl-pds">'</span></span>:</td> 1037 </tr> 1038 <tr> 1039 <td id="L155" class="blob-num js-line-number" data-line-number="155"></td> 1040 <td id="LC155" class="blob-code blob-code-inner js-file-line"> pars[<span class="pl-s"><span class="pl-pds">'</span>scale<span class="pl-pds">'</span></span>], pars[<span class="pl-s"><span class="pl-pds">'</span>background<span class="pl-pds">'</span></span>] <span class="pl-k">=</span> <span class="pl-c1">1</span>, <span class="pl-c1">0</span></td> 1041 </tr> 1042 <tr> 1043 <td id="L156" class="blob-num js-line-number" data-line-number="156"></td> 1044 <td id="LC156" class="blob-code blob-code-inner js-file-line"> 1045 </td> 1046 </tr> 1047 </table> 1048 1049 </div> 1050 1051 </div> 1052 1053 <a href="#jump-to-line" rel="facebox[.linejump]" data-hotkey="l" style="display:none">Jump to Line</a> 1054 <div id="jump-to-line" style="display:none"> 1055 <!-- </textarea> --><!-- '"` --><form accept-charset="UTF-8" action="" class="js-jump-to-line-form" method="get"><div style="margin:0;padding:0;display:inline"><input name="utf8" type="hidden" value="✓" /></div> 1056 <input class="linejump-input js-jump-to-line-field" type="text" placeholder="Jump to line…" aria-label="Jump to line" autofocus> 1057 <button type="submit" class="btn">Go</button> 1058 </form></div> 1059 1060 </div> 1061 <div class="modal-backdrop"></div> 1062 </div> 1063 1064 </div> 1065 </div> 1066 1067 </div> 1068 1069 <div class="container"> 1070 <div class="site-footer" role="contentinfo"> 1071 <ul class="site-footer-links right"> 1072 <li><a href="https://status.github.com/" data-ga-click="Footer, go to status, text:status">Status</a></li> 1073 <li><a href="https://developer.github.com" data-ga-click="Footer, go to api, text:api">API</a></li> 1074 <li><a href="https://training.github.com" data-ga-click="Footer, go to training, text:training">Training</a></li> 1075 <li><a href="https://shop.github.com" data-ga-click="Footer, go to shop, text:shop">Shop</a></li> 1076 <li><a href="https://github.com/blog" data-ga-click="Footer, go to blog, text:blog">Blog</a></li> 1077 <li><a href="https://github.com/about" data-ga-click="Footer, go to about, text:about">About</a></li> 1078 <li><a href="https://github.com/pricing" data-ga-click="Footer, go to pricing, text:pricing">Pricing</a></li> 1079 1080 </ul> 1081 1082 <a href="https://github.com" aria-label="Homepage"> 1083 <svg aria-hidden="true" class="octicon octicon-mark-github" height="24" role="img" title="GitHub " version="1.1" viewBox="0 0 16 16" width="24"><path d="M8 0C3.58 0 0 3.58 0 8c0 3.54 2.29 6.53 5.47 7.59 0.4 0.07 0.55-0.17 0.55-0.38 0-0.19-0.01-0.82-0.01-1.49-2.01 0.37-2.53-0.49-2.69-0.94-0.09-0.23-0.48-0.94-0.82-1.13-0.28-0.15-0.68-0.52-0.01-0.53 0.63-0.01 1.08 0.58 1.23 0.82 0.72 1.21 1.87 0.87 2.33 0.66 0.07-0.52 0.28-0.87 0.51-1.07-1.78-0.2-3.64-0.89-3.64-3.95 0-0.87 0.31-1.59 0.82-2.15-0.08-0.2-0.36-1.02 0.08-2.12 0 0 0.67-0.21 2.2 0.82 0.64-0.18 1.32-0.27 2-0.27 0.68 0 1.36 0.09 2 0.27 1.53-1.04 2.2-0.82 2.2-0.82 0.44 1.1 0.16 1.92 0.08 2.12 0.51 0.56 0.82 1.27 0.82 2.15 0 3.07-1.87 3.75-3.65 3.95 0.29 0.25 0.54 0.73 0.54 1.48 0 1.07-0.01 1.93-0.01 2.2 0 0.21 0.15 0.46 0.55 0.38C13.71 14.53 16 11.53 16 8 16 3.58 12.42 0 8 0z"></path></svg> 1084 </a> 1085 <ul class="site-footer-links"> 1086 <li>© 2016 <span title="0.05637s from github-fe127-cp1-prd.iad.github.net">GitHub</span>, Inc.</li> 1087 <li><a href="https://github.com/site/terms" data-ga-click="Footer, go to terms, text:terms">Terms</a></li> 1088 <li><a href="https://github.com/site/privacy" data-ga-click="Footer, go to privacy, text:privacy">Privacy</a></li> 1089 <li><a href="https://github.com/security" data-ga-click="Footer, go to security, text:security">Security</a></li> 1090 <li><a href="https://github.com/contact" data-ga-click="Footer, go to contact, text:contact">Contact</a></li> 1091 <li><a href="https://help.github.com" data-ga-click="Footer, go to help, text:help">Help</a></li> 1092 </ul> 1093 </div> 1094 </div> 1095 1096 1097 1098 1099 1100 1101 1102 <div id="ajax-error-message" class="flash flash-error"> 1103 <svg aria-hidden="true" class="octicon octicon-alert" height="16" role="img" version="1.1" viewBox="0 0 16 16" width="16"><path d="M15.72 12.5l-6.85-11.98C8.69 0.21 8.36 0.02 8 0.02s-0.69 0.19-0.87 0.5l-6.85 11.98c-0.18 0.31-0.18 0.69 0 1C0.47 13.81 0.8 14 1.15 14h13.7c0.36 0 0.69-0.19 0.86-0.5S15.89 12.81 15.72 12.5zM9 12H7V10h2V12zM9 9H7V5h2V9z"></path></svg> 1104 <button type="button" class="flash-close js-flash-close js-ajax-error-dismiss" aria-label="Dismiss error"> 1105 <svg aria-hidden="true" class="octicon octicon-x" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M7.48 8l3.75 3.75-1.48 1.48-3.75-3.75-3.75 3.75-1.48-1.48 3.75-3.75L0.77 4.25l1.48-1.48 3.75 3.75 3.75-3.75 1.48 1.48-3.75 3.75z"></path></svg> 1106 </button> 1107 Something went wrong with that request. Please try again. 1108 </div> 1109 1110 1111 1112 <script crossorigin="anonymous" integrity="sha256-7lIbjp+srGj/J+k/w64PjtgR17+eQ06E9LnqIneAsIQ=" src="https://assets-cdn.github.com/assets/frameworks-ee521b8e9facac68ff27e93fc3ae0f8ed811d7bf9e434e84f4b9ea227780b084.js"></script> 1113 <script async="async" crossorigin="anonymous" integrity="sha256-aWM2lkt8QuTG9N+68fjlf0Jc2anRihLz/m4910T9fRM=" src="https://assets-cdn.github.com/assets/github-696336964b7c42e4c6f4dfbaf1f8e57f425cd9a9d18a12f3fe6e3dd744fd7d13.js"></script> 1114 1115 1116 1117 <div class="js-stale-session-flash stale-session-flash flash flash-warn flash-banner hidden"> 1118 <svg aria-hidden="true" class="octicon octicon-alert" height="16" role="img" version="1.1" viewBox="0 0 16 16" width="16"><path d="M15.72 12.5l-6.85-11.98C8.69 0.21 8.36 0.02 8 0.02s-0.69 0.19-0.87 0.5l-6.85 11.98c-0.18 0.31-0.18 0.69 0 1C0.47 13.81 0.8 14 1.15 14h13.7c0.36 0 0.69-0.19 0.86-0.5S15.89 12.81 15.72 12.5zM9 12H7V10h2V12zM9 9H7V5h2V9z"></path></svg> 1119 <span class="signed-in-tab-flash">You signed in with another tab or window. <a href="">Reload</a> to refresh your session.</span> 1120 <span class="signed-out-tab-flash">You signed out in another tab or window. <a href="">Reload</a> to refresh your session.</span> 1121 </div> 1122 <div class="facebox" id="facebox" style="display:none;"> 1123 <div class="facebox-popup"> 1124 <div class="facebox-content" role="dialog" aria-labelledby="facebox-header" aria-describedby="facebox-description"> 1125 </div> 1126 <button type="button" class="facebox-close js-facebox-close" aria-label="Close modal"> 1127 <svg aria-hidden="true" class="octicon octicon-x" height="16" role="img" version="1.1" viewBox="0 0 12 16" width="12"><path d="M7.48 8l3.75 3.75-1.48 1.48-3.75-3.75-3.75 3.75-1.48-1.48 3.75-3.75L0.77 4.25l1.48-1.48 3.75 3.75 3.75-3.75 1.48 1.48-3.75 3.75z"></path></svg> 1128 </button> 1129 </div> 1130 </div> 1131 1132 </body> 1133 </html> 1134 -
sasmodels/core.py
rcde11f0f rd18582e 2 2 Core model handling routines. 3 3 """ 4 __all__ = ["list_models", "load_model_definition", "precompile_dll",5 "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR" ]6 4 7 5 from os.path import basename, dirname, join as joinpath … … 22 20 HAVE_OPENCL = False 23 21 22 __all__ = [ 23 "list_models", "load_model_definition", "precompile_dll", 24 "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR", 25 ] 24 26 25 27 def list_models(): … … 36 38 """ 37 39 Load a model definition given the model name. 40 41 This returns a handle to the module defining the model. This can be 42 used with functions in generate to build the docs or extract model info. 38 43 """ 39 44 __import__('sasmodels.models.'+model_name) … … 61 66 62 67 def isstr(s): 68 """ 69 Return True if *s* is a string-like object. 70 """ 63 71 try: s + '' 64 72 except: return False 65 73 return True 66 74 67 def load_model(model_definition, dtype= "single", platform="ocl"):75 def load_model(model_definition, dtype=None, platform="ocl"): 68 76 """ 69 77 Prepare the model for the default execution platform. … … 79 87 for the calculation. Any valid numpy single or double precision identifier 80 88 is valid, such as 'single', 'f', 'f32', or np.float32 for single, or 81 'double', 'd', 'f64' and np.float64 for double. 89 'double', 'd', 'f64' and np.float64 for double. If *None*, then use 90 'single' unless the model defines single=False. 82 91 83 92 *platform* should be "dll" to force the dll to be used for C models, … … 86 95 if isstr(model_definition): 87 96 model_definition = load_model_definition(model_definition) 97 if dtype is None: 98 dtype = 'single' if getattr(model_definition, 'single', True) else 'double' 88 99 source, info = generate.make(model_definition) 89 100 if callable(info.get('Iq', None)): … … 99 110 # source = open(info['name']+'.cl','r').read() 100 111 101 if (platform =="dll"112 if (platform == "dll" 102 113 or not HAVE_OPENCL 103 114 or not kernelcl.environment().has_type(dtype)): … … 110 121 Return a computation kernel from the model definition and the q input. 111 122 """ 112 model_input = model.make_input(q_vectors) 113 return model(model_input) 123 return model(q_vectors) 114 124 115 125 def get_weights(info, pars, name): … … 128 138 width = pars.get(name+'_pd', 0.0) 129 139 nsigma = pars.get(name+'_pd_nsigma', 3.0) 130 value, weight = weights.get_weights(140 value, weight = weights.get_weights( 131 141 disperser, npts, width, nsigma, value, limits, relative) 132 142 return value, weight / np.sum(weight) … … 195 205 for name in info['partype']['volume']] 196 206 value, weight = dispersion_mesh(vol_pars) 197 whole, part = VR(*value)207 whole, part = VR(*value) 198 208 return np.sum(weight*part)/np.sum(weight*whole) 199 209 -
sasmodels/data.py
r013adb7 rd18582e 87 87 88 88 class Data1D(object): 89 """ 90 1D data object. 91 92 Note that this definition matches the attributes from sasview, with 93 some generic 1D data vectors and some SAS specific definitions. Some 94 refactoring to allow consistent naming conventions between 1D, 2D and 95 SESANS data would be helpful. 96 97 **Attributes** 98 99 *x*, *dx*: $q$ vector and gaussian resolution 100 101 *y*, *dy*: $I(q)$ vector and measurement uncertainty 102 103 *mask*: values to include in plotting/analysis 104 105 *dxl*: slit widths for slit smeared data, with *dx* ignored 106 107 *qmin*, *qmax*: range of $q$ values in *x* 108 109 *filename*: label for the data line 110 111 *_xaxis*, *_xunit*: label and units for the *x* axis 112 113 *_yaxis*, *_yunit*: label and units for the *y* axis 114 """ 89 115 def __init__(self, x=None, y=None, dx=None, dy=None): 90 116 self.x, self.y, self.dx, self.dy = x, y, dx, dy 91 117 self.dxl = None 118 self.filename = None 119 self.qmin = x.min() if x is not None else np.NaN 120 self.qmax = x.max() if x is not None else np.NaN 121 # TODO: why is 1D mask False and 2D mask True? 122 self.mask = (np.isnan(y) if y is not None 123 else np.zeros_like(x, 'b') if x is not None 124 else None) 125 self._xaxis, self._xunit = "x", "" 126 self._yaxis, self._yunit = "y", "" 92 127 93 128 def xaxis(self, label, unit): … … 108 143 109 144 class Data2D(object): 110 def __init__(self): 145 """ 146 2D data object. 147 148 Note that this definition matches the attributes from sasview. Some 149 refactoring to allow consistent naming conventions between 1D, 2D and 150 SESANS data would be helpful. 151 152 **Attributes** 153 154 *qx_data*, *dqx_data*: $q_x$ matrix and gaussian resolution 155 156 *qy_data*, *dqy_data*: $q_y$ matrix and gaussian resolution 157 158 *data*, *err_data*: $I(q)$ matrix and measurement uncertainty 159 160 *mask*: values to exclude from plotting/analysis 161 162 *qmin*, *qmax*: range of $q$ values in *x* 163 164 *filename*: label for the data line 165 166 *_xaxis*, *_xunit*: label and units for the *x* axis 167 168 *_yaxis*, *_yunit*: label and units for the *y* axis 169 170 *_zaxis*, *_zunit*: label and units for the *y* axis 171 172 *Q_unit*, *I_unit*: units for Q and intensity 173 174 *x_bins*, *y_bins*: grid steps in *x* and *y* directions 175 """ 176 def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None): 177 self.qx_data, self.dqx_data = x, dx 178 self.qy_data, self.dqy_data = y, dy 179 self.data, self.err_data = z, dz 180 self.mask = (~np.isnan(z) if z is not None 181 else np.ones_like(x) if x is not None 182 else None) 183 self.q_data = np.sqrt(x**2 + y**2) 184 self.qmin = 1e-16 185 self.qmax = np.inf 111 186 self.detector = [] 112 187 self.source = Source() 188 self.Q_unit = "1/A" 189 self.I_unit = "1/cm" 190 self.xaxis("Q_x", "1/A") 191 self.yaxis("Q_y", "1/A") 192 self.zaxis("Intensity", "1/cm") 193 self._xaxis, self._xunit = "x", "" 194 self._yaxis, self._yunit = "y", "" 195 self._zaxis, self._zunit = "z", "" 196 self.x_bins, self.y_bins = None, None 113 197 114 198 def xaxis(self, label, unit): … … 135 219 136 220 class Vector(object): 221 """ 222 3-space vector of *x*, *y*, *z* 223 """ 137 224 def __init__(self, x=None, y=None, z=None): 138 225 self.x, self.y, self.z = x, y, z 139 226 140 227 class Detector(object): 228 """ 229 Detector attributes. 230 """ 231 def __init__(self, pixel_size=(None, None), distance=None): 232 self.pixel_size = Vector(*pixel_size) 233 self.distance = distance 234 235 class Source(object): 236 """ 237 Beam attributes. 238 """ 141 239 def __init__(self): 142 self.pixel_size = Vector() 143 144 class Source(object): 145 pass 146 147 148 def empty_data1D(q, resolution=0.05): 240 self.wavelength = np.NaN 241 self.wavelength_unit = "A" 242 243 244 def empty_data1D(q, resolution=0.0): 149 245 """ 150 246 Create empty 1D data using the given *q* as the x value. … … 156 252 #dIq = np.sqrt(Iq) 157 253 Iq, dIq = None, None 254 q = np.asarray(q) 158 255 data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 159 256 data.filename = "fake data" 160 data.qmin, data.qmax = q.min(), q.max()161 data.mask = np.zeros(len(q), dtype='bool')162 257 return data 163 258 164 259 165 def empty_data2D(qx, qy=None, resolution=0.0 5):260 def empty_data2D(qx, qy=None, resolution=0.0): 166 261 """ 167 262 Create empty 2D data using the given mesh. … … 173 268 if qy is None: 174 269 qy = qx 270 qx, qy = np.asarray(qx), np.asarray(qy) 271 # 5% dQ/Q resolution 175 272 Qx, Qy = np.meshgrid(qx, qy) 176 273 Qx, Qy = Qx.flatten(), Qy.flatten() 177 274 Iq = 100 * np.ones_like(Qx) 178 275 dIq = np.sqrt(Iq) 179 mask = np.ones(len(Iq), dtype='bool')180 181 data = Data2D()182 data.filename = "fake data"183 data.qx_data = Qx184 data.qy_data = Qy185 data.data = Iq186 data.err_data = dIq187 data.mask = mask188 data.qmin = 1e-16189 data.qmax = np.inf190 191 # 5% dQ/Q resolution192 276 if resolution != 0: 193 277 # https://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_15.pdf … … 197 281 # radial (which instead it should be inverse). 198 282 Q = np.sqrt(Qx**2 + Qy**2) 199 d ata.dqx_data= resolution * Q200 d ata.dqy_data= resolution * Q283 dqx = resolution * Q 284 dqy = resolution * Q 201 285 else: 202 data.dqx_data = data.dqy_data = None 203 204 detector = Detector() 205 detector.pixel_size.x = 5 # mm 206 detector.pixel_size.y = 5 # mm 207 detector.distance = 4 # m 208 data.detector.append(detector) 286 dqx = dqy = None 287 288 data = Data2D(x=Qx, y=Qy, z=Iq, dx=dqx, dy=dqy, dz=dIq) 209 289 data.x_bins = qx 210 290 data.y_bins = qy 291 data.filename = "fake data" 292 293 # pixel_size in mm, distance in m 294 detector = Detector(pixel_size=(5, 5), distance=4) 295 data.detector.append(detector) 211 296 data.source.wavelength = 5 # angstroms 212 297 data.source.wavelength_unit = "A" 213 data.Q_unit = "1/A"214 data.I_unit = "1/cm"215 data.q_data = np.sqrt(Qx ** 2 + Qy ** 2)216 data.xaxis("Q_x", "A^{-1}")217 data.yaxis("Q_y", "A^{-1}")218 data.zaxis("Intensity", r"\text{cm}^{-1}")219 298 return data 220 299 … … 223 302 """ 224 303 Plot data loaded by the sasview loader. 304 305 *data* is a sasview data object, either 1D, 2D or SESANS. 306 307 *view* is log or linear. 308 309 *limits* sets the intensity limits on the plot; if None then the limits 310 are inferred from the data. 225 311 """ 226 312 # Note: kind of weird using the plot result functions to plot just the … … 228 314 # do not repeat. 229 315 if hasattr(data, 'lam'): 230 _plot_result_sesans(data, None, None, plot_data=True, limits=limits)316 _plot_result_sesans(data, None, None, use_data=True, limits=limits) 231 317 elif hasattr(data, 'qx_data'): 232 _plot_result2D(data, None, None, view, plot_data=True, limits=limits)318 _plot_result2D(data, None, None, view, use_data=True, limits=limits) 233 319 else: 234 _plot_result1D(data, None, None, view, plot_data=True, limits=limits)320 _plot_result1D(data, None, None, view, use_data=True, limits=limits) 235 321 236 322 237 323 def plot_theory(data, theory, resid=None, view='log', 238 plot_data=True, limits=None): 324 use_data=True, limits=None): 325 """ 326 Plot theory calculation. 327 328 *data* is needed to define the graph properties such as labels and 329 units, and to define the data mask. 330 331 *theory* is a matrix of the same shape as the data. 332 333 *view* is log or linear 334 335 *use_data* is True if the data should be plotted as well as the theory. 336 337 *limits* sets the intensity limits on the plot; if None then the limits 338 are inferred from the data. 339 """ 239 340 if hasattr(data, 'lam'): 240 _plot_result_sesans(data, theory, resid, plot_data=True, limits=limits)341 _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 241 342 elif hasattr(data, 'qx_data'): 242 _plot_result2D(data, theory, resid, view, plot_data, limits=limits)343 _plot_result2D(data, theory, resid, view, use_data, limits=limits) 243 344 else: 244 _plot_result1D(data, theory, resid, view, plot_data, limits=limits)345 _plot_result1D(data, theory, resid, view, use_data, limits=limits) 245 346 246 347 247 348 def protect(fn): 349 """ 350 Decorator to wrap calls in an exception trapper which prints the 351 exception and continues. Keyboard interrupts are ignored. 352 """ 248 353 def wrapper(*args, **kw): 354 """ 355 Trap and print errors from function. 356 """ 249 357 try: 250 358 return fn(*args, **kw) 359 except KeyboardInterrupt: 360 raise 251 361 except: 252 362 traceback.print_exc() 253 pass254 363 255 364 return wrapper … … 257 366 258 367 @protect 259 def _plot_result1D(data, theory, resid, view, plot_data, limits=None):368 def _plot_result1D(data, theory, resid, view, use_data, limits=None): 260 369 """ 261 370 Plot the data and residuals for 1D data. … … 264 373 from numpy.ma import masked_array, masked 265 374 266 plot_theory = theory is not None267 plot_resid = residis not None268 269 if data.y is None:270 plot_data = False 375 use_data = use_data and data.y is not None 376 use_theory = theory is not None 377 use_resid = resid is not None 378 num_plots = (use_data or use_theory) + use_resid 379 271 380 scale = data.x**4 if view == 'q4' else 1.0 272 381 273 if plot_data or plot_theory: 274 if plot_resid: 275 plt.subplot(121) 276 382 if use_data or use_theory: 277 383 #print(vmin, vmax) 278 384 all_positive = True 279 385 some_present = False 280 if plot_data:386 if use_data: 281 387 mdata = masked_array(data.y, data.mask.copy()) 282 388 mdata[~np.isfinite(mdata)] = masked … … 284 390 mdata[mdata <= 0] = masked 285 391 plt.errorbar(data.x/10, scale*mdata, yerr=data.dy, fmt='.') 286 all_positive = all_positive and (mdata >0).all()392 all_positive = all_positive and (mdata > 0).all() 287 393 some_present = some_present or (mdata.count() > 0) 288 394 289 395 290 if plot_theory:396 if use_theory: 291 397 mtheory = masked_array(theory, data.mask.copy()) 292 398 mtheory[~np.isfinite(mtheory)] = masked 293 399 if view is 'log': 294 mtheory[mtheory <=0] = masked400 mtheory[mtheory <= 0] = masked 295 401 plt.plot(data.x/10, scale*mtheory, '-', hold=True) 296 all_positive = all_positive and (mtheory >0).all()402 all_positive = all_positive and (mtheory > 0).all() 297 403 some_present = some_present or (mtheory.count() > 0) 298 404 299 405 if limits is not None: 300 406 plt.ylim(*limits) 407 408 if num_plots > 1: 409 plt.subplot(1, num_plots, 1) 301 410 plt.xscale('linear' if not some_present else view) 302 411 plt.yscale('linear' … … 306 415 plt.ylabel('$I(q)$') 307 416 308 if plot_resid: 309 if plot_data or plot_theory: 310 plt.subplot(122) 311 417 if use_resid: 312 418 mresid = masked_array(resid, data.mask.copy()) 313 419 mresid[~np.isfinite(mresid)] = masked 314 420 some_present = (mresid.count() > 0) 421 422 if num_plots > 1: 423 plt.subplot(1, num_plots, (use_data or use_theory) + 1) 315 424 plt.plot(data.x/10, mresid, '-') 316 425 plt.xlabel("$q$/nm$^{-1}$") … … 320 429 321 430 @protect 322 def _plot_result_sesans(data, theory, resid, plot_data, limits=None): 431 def _plot_result_sesans(data, theory, resid, use_data, limits=None): 432 """ 433 Plot SESANS results. 434 """ 323 435 import matplotlib.pyplot as plt 324 if data.y is None:325 plot_data = False326 plot_theory = theoryis not None327 plot_resid = resid is not None328 329 if plot_data or plot_theory:330 if plot_resid:331 plt.subplot(1 21)332 if plot_data:436 use_data = use_data and data.y is not None 437 use_theory = theory is not None 438 use_resid = resid is not None 439 num_plots = (use_data or use_theory) + use_resid 440 441 if use_data or use_theory: 442 if num_plots > 1: 443 plt.subplot(1, num_plots, 1) 444 if use_data: 333 445 plt.errorbar(data.x, data.y, yerr=data.dy) 334 446 if theory is not None: … … 340 452 341 453 if resid is not None: 342 if plot_data or plot_theory: 343 plt.subplot(122) 344 454 if num_plots > 1: 455 plt.subplot(1, num_plots, (use_data or use_theory) + 1) 345 456 plt.plot(data.x, resid, 'x') 346 457 plt.xlabel('spin echo length (nm)') … … 349 460 350 461 @protect 351 def _plot_result2D(data, theory, resid, view, plot_data, limits=None):462 def _plot_result2D(data, theory, resid, view, use_data, limits=None): 352 463 """ 353 464 Plot the data and residuals for 2D data. 354 465 """ 355 466 import matplotlib.pyplot as plt 356 if data.data is None:357 plot_data = False358 plot_theory = theoryis not None359 plot_resid = resid is not None467 use_data = use_data and data.data is not None 468 use_theory = theory is not None 469 use_resid = resid is not None 470 num_plots = use_data + use_theory + use_resid 360 471 361 472 # Put theory and data on a common colormap scale 362 if limits is None: 363 vmin, vmax = np.inf, -np.inf 364 if plot_data: 365 target = data.data[~data.mask] 366 datamin = target[target>0].min() if view == 'log' else target.min() 367 datamax = target.max() 368 vmin = min(vmin, datamin) 369 vmax = max(vmax, datamax) 370 if plot_theory: 371 theorymin = theory[theory>0].min() if view=='log' else theory.min() 372 theorymax = theory.max() 373 vmin = min(vmin, theorymin) 374 vmax = max(vmax, theorymax) 375 else: 473 vmin, vmax = np.inf, -np.inf 474 if use_data: 475 target = data.data[~data.mask] 476 datamin = target[target > 0].min() if view == 'log' else target.min() 477 datamax = target.max() 478 vmin = min(vmin, datamin) 479 vmax = max(vmax, datamax) 480 if use_theory: 481 theorymin = theory[theory > 0].min() if view == 'log' else theory.min() 482 theorymax = theory.max() 483 vmin = min(vmin, theorymin) 484 vmax = max(vmax, theorymax) 485 486 # Override data limits from the caller 487 if limits is not None: 376 488 vmin, vmax = limits 377 489 378 if plot_data: 379 if plot_theory and plot_resid: 380 plt.subplot(131) 381 elif plot_theory or plot_resid: 382 plt.subplot(121) 490 # Plot data 491 if use_data: 492 if num_plots > 1: 493 plt.subplot(1, num_plots, 1) 383 494 _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax) 384 495 plt.title('data') … … 386 497 h.set_label('$I(q)$') 387 498 388 if plot_theory: 389 if plot_data and plot_resid: 390 plt.subplot(132) 391 elif plot_data: 392 plt.subplot(122) 393 elif plot_resid: 394 plt.subplot(121) 499 # plot theory 500 if use_theory: 501 if num_plots > 1: 502 plt.subplot(1, num_plots, use_data+1) 395 503 _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax) 396 504 plt.title('theory') 397 505 h = plt.colorbar() 398 h.set_label(r'$\log_{10}I(q)$' if view =='log'506 h.set_label(r'$\log_{10}I(q)$' if view == 'log' 399 507 else r'$q^4 I(q)$' if view == 'q4' 400 508 else '$I(q)$') 401 509 402 #if plot_data or plot_theory: 403 # plt.colorbar() 404 405 if plot_resid: 406 if plot_data and plot_theory: 407 plt.subplot(133) 408 elif plot_data or plot_theory: 409 plt.subplot(122) 510 # plot resid 511 if use_resid: 512 if num_plots > 1: 513 plt.subplot(1, num_plots, use_data+use_theory+1) 410 514 _plot_2d_signal(data, resid, view='linear') 411 515 plt.title('residuals') 412 516 h = plt.colorbar() 413 h.set_label( '$\Delta I(q)$')517 h.set_label(r'$\Delta I(q)$') 414 518 415 519 … … 456 560 457 561 def demo(): 562 """ 563 Load and plot a SAS dataset. 564 """ 458 565 data = load_data('DEC07086.DAT') 459 566 set_beam_stop(data, 0.004) -
sasmodels/direct_model.py
r9404dd3 rd18582e 1 import warnings 1 """ 2 Class interface to the model calculator. 3 4 Calling a model is somewhat non-trivial since the functions called depend 5 on the data type. For 1D data the *Iq* kernel needs to be called, for 6 2D data the *Iqxy* kernel needs to be called, and for SESANS data the 7 *Iq* kernel needs to be called followed by a Hankel transform. Before 8 the kernel is called an appropriate *q* calculation vector needs to be 9 constructed. This is not the simple *q* vector where you have measured 10 the data since the resolution calculation will require values beyond the 11 range of the measured data. After the calculation the resolution calculator 12 must be called to return the predicted value for each measured data point. 13 14 :class:`DirectModel` is a callable object that takes *parameter=value* 15 keyword arguments and returns the appropriate theory values for the data. 16 17 :class:`DataMixin` does the real work of interpreting the data and calling 18 the model calculator. This is used by :class:`DirectModel`, which uses 19 direct parameter values and by :class:`bumps_model.Experiment` which wraps 20 the parameter values in boxes so that the user can set fitting ranges, etc. 21 on the individual parameters and send the model to the Bumps optimizers. 22 """ 23 from __future__ import print_function 2 24 3 25 import numpy as np … … 19 41 currently used by :class:`sasview_model.SasviewModel` since this will 20 42 require a number of changes to SasView before we can do it. 43 44 :meth:`_interpret_data` initializes the data structures necessary 45 to manage the calculations. This sets attributes in the child class 46 such as *data_type* and *resolution*. 47 48 :meth:`_calc_theory` evaluates the model at the given control values. 49 50 :meth:`_set_data` sets the intensity data in the data object, 51 possibly with random noise added. This is useful for simulating a 52 dataset with the results from :meth:`_calc_theory`. 21 53 """ 22 54 def _interpret_data(self, data, model): 55 # pylint: disable=attribute-defined-outside-init 56 23 57 self._data = data 24 58 self._model = model … … 36 70 if self.data_type == 'sesans': 37 71 q = sesans.make_q(data.sample.zacceptance, data.Rmax) 38 self.index = slice(None, None) 72 index = slice(None, None) 73 res = None 39 74 if data.y is not None: 40 self.Iq = data.y 41 self.dIq = data.dy 75 Iq, dIq = data.y, data.dy 76 else: 77 Iq, dIq = None, None 42 78 #self._theory = np.zeros_like(q) 43 79 q_vectors = [q] … … 49 85 qmax = getattr(data, 'qmax', np.inf) 50 86 accuracy = getattr(data, 'accuracy', 'Low') 51 self.index = ~data.mask & (q >= qmin) & (q <= qmax)87 index = ~data.mask & (q >= qmin) & (q <= qmax) 52 88 if data.data is not None: 53 self.index &= ~np.isnan(data.data) 54 self.Iq = data.data[self.index] 55 self.dIq = data.err_data[self.index] 56 self.resolution = resolution2d.Pinhole2D(data=data, index=self.index, 57 nsigma=3.0, accuracy=accuracy) 89 index &= ~np.isnan(data.data) 90 Iq = data.data[index] 91 dIq = data.err_data[index] 92 else: 93 Iq, dIq = None, None 94 res = resolution2d.Pinhole2D(data=data, index=index, 95 nsigma=3.0, accuracy=accuracy) 58 96 #self._theory = np.zeros_like(self.Iq) 59 q_vectors = self.resolution.q_calc97 q_vectors = res.q_calc 60 98 elif self.data_type == 'Iq': 61 self.index = (data.x >= data.qmin) & (data.x <= data.qmax)99 index = (data.x >= data.qmin) & (data.x <= data.qmax) 62 100 if data.y is not None: 63 self.index &= ~np.isnan(data.y) 64 self.Iq = data.y[self.index] 65 self.dIq = data.dy[self.index] 101 index &= ~np.isnan(data.y) 102 Iq = data.y[index] 103 dIq = data.dy[index] 104 else: 105 Iq, dIq = None, None 66 106 if getattr(data, 'dx', None) is not None: 67 q, dq = data.x[ self.index], data.dx[self.index]68 if (dq >0).any():69 self.resolution= resolution.Pinhole1D(q, dq)107 q, dq = data.x[index], data.dx[index] 108 if (dq > 0).any(): 109 res = resolution.Pinhole1D(q, dq) 70 110 else: 71 self.resolution= resolution.Perfect1D(q)72 elif (getattr(data, 'dxl', None) is not None and73 74 self.resolution = resolution.Slit1D(data.x[self.index],75 width=data.dxh[self.index],76 height=data.dxw[self.index])77 else: 78 self.resolution = resolution.Perfect1D(data.x[self.index])111 res = resolution.Perfect1D(q) 112 elif (getattr(data, 'dxl', None) is not None 113 and getattr(data, 'dxw', None) is not None): 114 res = resolution.Slit1D(data.x[index], 115 width=data.dxh[index], 116 height=data.dxw[index]) 117 else: 118 res = resolution.Perfect1D(data.x[index]) 79 119 80 120 #self._theory = np.zeros_like(self.Iq) 81 q_vectors = [ self.resolution.q_calc]121 q_vectors = [res.q_calc] 82 122 else: 83 123 raise ValueError("Unknown data type") # never gets here … … 87 127 self._kernel_inputs = [v for v in q_vectors] 88 128 self._kernel = None 129 self.Iq, self.dIq, self.index = Iq, dIq, index 130 self.resolution = res 89 131 90 132 def _set_data(self, Iq, noise=None): 133 # pylint: disable=attribute-defined-outside-init 91 134 if noise is not None: 92 135 self.dIq = Iq*noise*0.01 … … 106 149 def _calc_theory(self, pars, cutoff=0.0): 107 150 if self._kernel is None: 108 q_input = self._model.make_input(self._kernel_inputs) 109 self._kernel = self._model(q_input) 151 self._kernel = make_kernel(self._model, self._kernel_inputs) # pylint: disable=attribute-defined-outside-init 110 152 111 153 Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) … … 120 162 121 163 class DirectModel(DataMixin): 164 """ 165 Create a calculator object for a model. 166 167 *data* is 1D SAS, 2D SAS or SESANS data 168 169 *model* is a model calculator return from :func:`generate.load_model` 170 171 *cutoff* is the polydispersity weight cutoff. 172 """ 122 173 def __init__(self, data, model, cutoff=1e-5): 123 174 self.model = model 124 175 self.cutoff = cutoff 176 # Note: _interpret_data defines the model attributes 125 177 self._interpret_data(data, model) 126 self.kernel = make_kernel(self.model, self._kernel_inputs) 178 127 179 def __call__(self, **pars): 128 180 return self._calc_theory(pars, cutoff=self.cutoff) 181 129 182 def ER(self, **pars): 183 """ 184 Compute the equivalent radius for the model. 185 186 Return 0. if not defined. 187 """ 130 188 return call_ER(self.model.info, pars) 189 131 190 def VR(self, **pars): 191 """ 192 Compute the equivalent volume for the model, including polydispersity 193 effects. 194 195 Return 1. if not defined. 196 """ 132 197 return call_VR(self.model.info, pars) 198 133 199 def simulate_data(self, noise=None, **pars): 200 """ 201 Generate simulated data for the model. 202 """ 134 203 Iq = self.__call__(**pars) 135 204 self._set_data(Iq, noise=noise) 136 205 137 def demo(): 206 def main(): 207 """ 208 Program to evaluate a particular model at a set of q values. 209 """ 138 210 import sys 139 211 from .data import empty_data1D, empty_data2D … … 144 216 model_name = sys.argv[1] 145 217 call = sys.argv[2].upper() 146 if call not in ("ER", "VR"):218 if call not in ("ER", "VR"): 147 219 try: 148 220 values = [float(v) for v in call.split(',')] … … 153 225 data = empty_data1D([q]) 154 226 elif len(values) == 2: 155 qx, qy = values156 data = empty_data2D([qx], [qy])227 qx, qy = values 228 data = empty_data2D([qx], [qy]) 157 229 else: 158 230 print("use q or qx,qy or ER or VR") … … 162 234 163 235 model_definition = load_model_definition(model_name) 164 model = load_model(model_definition , dtype='single')236 model = load_model(model_definition) 165 237 calculator = DirectModel(data, model) 166 pars = dict((k, float(v))238 pars = dict((k, float(v)) 167 239 for pair in sys.argv[3:] 168 for k, v in [pair.split('=')])240 for k, v in [pair.split('=')]) 169 241 if call == "ER": 170 242 print(calculator.ER(**pars)) … … 176 248 177 249 if __name__ == "__main__": 178 demo()250 main() -
sasmodels/exception.py
reaca9eb r823e620 8 8 except NameError: 9 9 class WindowsError(Exception): 10 """ 11 Fake WindowsException when not on Windows. 12 """ 10 13 pass 11 14 … … 31 34 # TODO: try to incorporate current stack trace in the raised exception 32 35 if isinstance(exc, WindowsError): 33 raise OSError(str(exc) +" "+msg)36 raise OSError(str(exc) + " " + msg) 34 37 35 38 args = exc.args … … 38 41 else: 39 42 try: 40 arg0 = " ".join((args[0], msg))43 arg0 = " ".join((args[0], msg)) 41 44 exc.args = tuple([arg0] + list(args[1:])) 42 45 except: 43 exc.args = (" ".join((str(exc), msg)),)46 exc.args = (" ".join((str(exc), msg)),) -
sasmodels/generate.py
re66c9f9 reafc9fa 7 7 particular dimensions averaged over all orientations. 8 8 9 *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form9 *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 10 10 with particular dimensions for a single orientation. 11 11 … … 47 47 the *sin* and *cos* values are needed for a particular argument. Since 48 48 this function does not exist in C99, all use of *sincos* should be 49 replaced by the macro *SINCOS(value, sn,cn)* where *sn* and *cn* are49 replaced by the macro *SINCOS(value, sn, cn)* where *sn* and *cn* are 50 50 previously declared *double* variables. When compiled for systems without 51 51 OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls. If *value* is … … 194 194 returns a list of files required by the model. 195 195 """ 196 from __future__ import print_function 196 197 197 198 # TODO: identify model files which have changed since loading and reload them. 198 199 __all__ = ["make", "doc", "sources", "convert_type"]200 199 201 200 import sys … … 206 205 207 206 import numpy as np 207 208 __all__ = ["make", "doc", "sources", "convert_type"] 209 208 210 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 209 211 … … 216 218 F128 = None 217 219 218 219 220 # Scale and background, which are parameters common to every form factor 220 221 COMMON_PARAMETERS = [ … … 222 223 ["background", "1/cm", 0, [0, np.inf], "", "Source background"], 223 224 ] 224 225 225 226 226 # Conversion from units defined in the parameter table for each model … … 266 266 267 267 def format_units(units): 268 """ 269 Convert units into ReStructured Text format. 270 """ 268 271 return "string" if isinstance(units, list) else RST_UNITS.get(units, units) 269 272 … … 310 313 raise ValueError("%r not found in %s" % (filename, search_path)) 311 314 312 def sources(info):315 def model_sources(info): 313 316 """ 314 317 Return a list of the sources file paths for the module. … … 335 338 """ 336 339 Convert code from double precision to the desired type. 340 341 Floating point constants are tagged with 'f' for single precision or 'L' 342 for long double precision. 337 343 """ 338 344 if dtype == F16: … … 350 356 351 357 def _convert_type(source, type_name, constant_flag): 358 """ 359 Replace 'double' with *type_name* in *source*, tagging floating point 360 constants with *constant_flag*. 361 """ 352 362 # Convert double keyword to float/long double/half. 353 363 # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16. … … 363 373 364 374 365 def kernel_name(info, is_2 D):375 def kernel_name(info, is_2d): 366 376 """ 367 377 Name of the exported kernel symbol. 368 378 """ 369 return info['name'] + "_" + ("Iqxy" if is_2 Delse "Iq")379 return info['name'] + "_" + ("Iqxy" if is_2d else "Iq") 370 380 371 381 … … 410 420 411 421 412 def build_polydispersity_loops(pd_pars): 413 """ 414 Build polydispersity loops 415 416 Returns loop opening and loop closing 417 """ 418 LOOP_OPEN = """\ 422 LOOP_OPEN = """\ 419 423 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 420 424 const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 421 425 const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 422 426 """ 427 def build_polydispersity_loops(pd_pars): 428 """ 429 Build polydispersity loops 430 431 Returns loop opening and loop closing 432 """ 423 433 depth = 4 424 434 offset = "" … … 455 465 456 466 # Load additional sources 457 source = [open(f).read() for f in sources(info)]467 source = [open(f).read() for f in model_sources(info)] 458 468 459 469 # Prepare defines … … 563 573 564 574 #for d in defines: print(d) 565 DEFINES= '\n'.join('#define %s %s' % (k, v) for k, v in defines)566 SOURCES= '\n\n'.join(source)575 defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 576 sources = '\n\n'.join(source) 567 577 return C_KERNEL_TEMPLATE % { 568 'DEFINES': DEFINES,569 'SOURCES': SOURCES,578 'DEFINES': defines, 579 'SOURCES': sources, 570 580 } 571 581 … … 581 591 demo_parameters = getattr(kernel_module, 'demo', None) 582 592 if demo_parameters is None: 583 demo_parameters = dict((p[0], p[2]) for p in parameters)593 demo_parameters = dict((p[0], p[2]) for p in parameters) 584 594 filename = abspath(kernel_module.__file__) 585 595 kernel_id = splitext(basename(filename))[0] … … 588 598 name = " ".join(w.capitalize() for w in kernel_id.split('_')) 589 599 info = dict( 590 id =kernel_id, # string used to load the kernel600 id=kernel_id, # string used to load the kernel 591 601 filename=abspath(kernel_module.__file__), 592 602 name=name, … … 625 635 %re.escape(string.punctuation)) 626 636 def _convert_section_titles_to_boldface(lines): 637 """ 638 Do the actual work of identifying and converting section headings. 639 """ 627 640 prior = None 628 641 for line in lines: … … 631 644 elif section_marker.match(line): 632 645 if len(line) >= len(prior): 633 yield "".join( ("**",prior,"**"))646 yield "".join(("**", prior, "**")) 634 647 prior = None 635 648 else: … … 642 655 yield prior 643 656 644 def convert_section_titles_to_boldface(string): 645 return "\n".join(_convert_section_titles_to_boldface(string.split('\n'))) 657 def convert_section_titles_to_boldface(s): 658 """ 659 Use explicit bold-face rather than section headings so that the table of 660 contents is not polluted with section names from the model documentation. 661 662 Sections are identified as the title line followed by a line of punctuation 663 at least as long as the title line. 664 """ 665 return "\n".join(_convert_section_titles_to_boldface(s.split('\n'))) 646 666 647 667 def doc(kernel_module): … … 666 686 667 687 def demo_time(): 688 """ 689 Show how long it takes to process a model. 690 """ 668 691 from .models import cylinder 669 692 import datetime … … 674 697 675 698 def main(): 699 """ 700 Program which prints the source produced by the model. 701 """ 676 702 if len(sys.argv) <= 1: 677 703 print("usage: python -m sasmodels.generate modelname") -
sasmodels/kernel_template.c
r9c79c32 r840b859 16 16 using namespace std; 17 17 #if defined(_MSC_VER) 18 # define kernel extern "C" __declspec( dllexport ) 18 #include <limits> 19 #include <float.h> 20 #define kernel extern "C" __declspec( dllexport ) 19 21 inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } 20 inline double fmin(double x, double y) { return x>y ? y : x; } 21 inline double fmax(double x, double y) { return x<y ? y : x; } 22 inline double fmin(double x, double y) { return x>y ? y : x; } 23 inline double fmax(double x, double y) { return x<y ? y : x; } 24 inline double isnan(double x) { return _isnan(x); } 25 #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN 26 static double cephes_expm1(double x) { 27 // Adapted from the cephes math library. 28 // Copyright 1984 - 1992 by Stephen L. Moshier 29 if (x != x || x == 0.0) { 30 return x; // NaN and +/- 0 31 } else if (x < -0.5 || x > 0.5) { 32 return exp(x) - 1.0; 33 } else { 34 const double xsq = x*x; 35 const double p = ((( 36 +1.2617719307481059087798E-4)*xsq 37 +3.0299440770744196129956E-2)*xsq 38 +9.9999999999999999991025E-1); 39 const double q = (((( 40 +3.0019850513866445504159E-6)*xsq 41 +2.5244834034968410419224E-3)*xsq 42 +2.2726554820815502876593E-1)*xsq 43 +2.0000000000000000000897E0); 44 double r = x * p; 45 r = r / (q - r); 46 return r+r; 47 } 48 } 49 #define expm1 cephes_expm1 22 50 #else 23 #define kernel extern "C"51 #define kernel extern "C" 24 52 #endif 25 53 inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); } -
sasmodels/kernelcl.py
rcde11f0f re6a5556 1 1 """ 2 GPU support through OpenCL2 GPU driver for C kernels 3 3 4 4 There should be a single GPU environment running on the system. This … … 152 152 153 153 154 def make_result(self, size):155 self.res = np.empty(size, dtype=self.dtype)156 self.res_b = cl.Buffer(self.program.context, mf.READ_WRITE, self.res.nbytes)157 return self.res, self.res_b158 159 160 154 # for now, this returns one device in the context 161 155 # TODO: create a context that contains all devices on all platforms … … 178 172 #self.data_boundary = max(d.min_data_type_align_size 179 173 # for d in self.context.devices) 180 self.queues = [cl.CommandQueue( self.context, d)181 for d in self.context.devices]174 self.queues = [cl.CommandQueue(context, context.devices[0]) 175 for context in self.context] 182 176 self.compiled = {} 183 177 184 178 def has_type(self, dtype): 179 """ 180 Return True if all devices support a given type. 181 """ 185 182 dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype) 186 return all(has_type(d, dtype) for d in self.context.devices) 183 return any(has_type(d, dtype) 184 for context in self.context 185 for d in context.devices) 186 187 def get_queue(self, dtype): 188 """ 189 Return a command queue for the kernels of type dtype. 190 """ 191 for context, queue in zip(self.context, self.queues): 192 if all(has_type(d, dtype) for d in context.devices): 193 return queue 194 195 def get_context(self, dtype): 196 """ 197 Return a OpenCL context for the kernels of type dtype. 198 """ 199 for context, queue in zip(self.context, self.queues): 200 if all(has_type(d, dtype) for d in context.devices): 201 return context 187 202 188 203 def _create_some_context(self): 204 """ 205 Protected call to cl.create_some_context without interactivity. Use 206 this if PYOPENCL_CTX is set in the environment. Sets the *context* 207 attribute. 208 """ 189 209 try: 190 self.context = cl.create_some_context(interactive=False)210 self.context = [cl.create_some_context(interactive=False)] 191 211 except Exception as exc: 192 212 warnings.warn(str(exc)) … … 195 215 196 216 def compile_program(self, name, source, dtype, fast=False): 217 """ 218 Compile the program for the device in the given context. 219 """ 197 220 key = "%s-%s-%s"%(name, dtype, fast) 198 221 if key not in self.compiled: 199 222 #print("compiling",name) 200 223 dtype = np.dtype(dtype) 201 program = compile_model(self. context, source, dtype, fast)224 program = compile_model(self.get_context(dtype), source, dtype, fast) 202 225 self.compiled[key] = program 203 226 return self.compiled[key] 204 227 205 228 def release_program(self, name): 229 """ 230 Free memory associated with the program on the device. 231 """ 206 232 if name in self.compiled: 207 233 self.compiled[name].release() … … 209 235 210 236 def _get_default_context(): 211 default = None 237 """ 238 Get an OpenCL context, preferring GPU over CPU, and preferring Intel 239 drivers over AMD drivers. 240 """ 241 # Note: on mobile devices there is automatic clock scaling if either the 242 # CPU or the GPU is underutilized; probably doesn't affect us, but we if 243 # it did, it would mean that putting a busy loop on the CPU while the GPU 244 # is running may increase throughput. 245 # 246 # Macbook pro, base install: 247 # {'Apple': [Intel CPU, NVIDIA GPU]} 248 # Macbook pro, base install: 249 # {'Apple': [Intel CPU, Intel GPU]} 250 # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed 251 # {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]} 252 gpu, cpu = None, None 212 253 for platform in cl.get_platforms(): 254 # AMD provides a much weaker CPU driver than Intel/Apple, so avoid it. 255 # If someone has bothered to install the AMD/NVIDIA drivers, prefer them over the integrated 256 # graphics driver that may have been supplied with the CPU chipset. 257 preferred_cpu = platform.vendor.startswith('Intel') or platform.vendor.startswith('Apple') 258 preferred_gpu = platform.vendor.startswith('Advanced') or platform.vendor.startswith('NVIDIA') 213 259 for device in platform.get_devices(): 214 260 if device.type == cl.device_type.GPU: 215 return cl.Context([device]) 216 if default is None: 217 default = device 218 219 if not default: 220 raise RuntimeError("OpenCL device not found") 221 222 return cl.Context([default]) 261 # If the existing type is not GPU then it will be CUSTOM or ACCELERATOR, 262 # so don't override it. 263 if gpu is None or (preferred_gpu and gpu.type == cl.device_type.GPU): 264 gpu = device 265 elif device.type == cl.device_type.CPU: 266 if cpu is None or preferred_cpu: 267 cpu = device 268 else: 269 # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM 270 # Intel Phi for example registers as an accelerator 271 # Since the user installed a custom device on their system and went through the 272 # pain of sorting out OpenCL drivers for it, lets assume they really do want to 273 # use it as their primary compute device. 274 gpu = device 275 276 # order the devices by gpu then by cpu; when searching for an available device by data type they 277 # will be checked in this order, which means that if the gpu supports double then the cpu will never 278 # be used (though we may make it possible to explicitly request the cpu at some point). 279 devices = [] 280 if gpu is not None: 281 devices.append(gpu) 282 if cpu is not None: 283 devices.append(cpu) 284 return [cl.Context([d]) for d in devices] 223 285 224 286 … … 241 303 self.info = info 242 304 self.source = source 243 self.dtype = generate.F32 if dtype =='fast' else np.dtype(dtype)305 self.dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype) 244 306 self.fast = (dtype == 'fast') 245 307 self.program = None # delay program creation 246 308 247 309 def __getstate__(self): 248 state = self.__dict__.copy() 249 state['program'] = None 250 return state 310 return self.info, self.source, self.dtype, self.fast 251 311 252 312 def __setstate__(self, state): 253 self.__dict__ = state.copy() 254 255 def __call__(self, q_input): 256 if self.dtype != q_input.dtype: 257 raise TypeError("data is %s kernel is %s" 258 % (q_input.dtype, self.dtype)) 313 self.info, self.source, self.dtype, self.fast = state 314 self.program = None 315 316 def __call__(self, q_vectors): 259 317 if self.program is None: 260 318 compiler = environment().compile_program 261 319 self.program = compiler(self.info['name'], self.source, self.dtype, 262 320 self.fast) 263 kernel_name = generate.kernel_name(self.info, q_input.is_2D) 321 is_2d = len(q_vectors) == 2 322 kernel_name = generate.kernel_name(self.info, is_2d) 264 323 kernel = getattr(self.program, kernel_name) 265 return GpuKernel(kernel, self.info, q_ input)324 return GpuKernel(kernel, self.info, q_vectors, self.dtype) 266 325 267 326 def release(self): 327 """ 328 Free the resources associated with the model. 329 """ 268 330 if self.program is not None: 269 331 environment().release_program(self.info['name']) 270 332 self.program = None 271 333 272 def make_input(self, q_vectors): 273 """ 274 Make q input vectors available to the model. 275 276 Note that each model needs its own q vector even if the case of 277 mixture models because some models may be OpenCL, some may be 278 ctypes and some may be pure python. 279 """ 280 return GpuInput(q_vectors, dtype=self.dtype) 334 def __del__(self): 335 self.release() 281 336 282 337 # TODO: check that we don't need a destructor for buffers which go out of scope … … 304 359 self.nq = q_vectors[0].size 305 360 self.dtype = np.dtype(dtype) 306 self.is_2 D= (len(q_vectors) == 2)361 self.is_2d = (len(q_vectors) == 2) 307 362 # TODO: stretch input based on get_warp() 308 363 # not doing it now since warp depends on kernel, which is not known … … 310 365 # architectures tested so far. 311 366 self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 367 context = env.get_context(self.dtype) 312 368 self.q_buffers = [ 313 cl.Buffer( env.context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)369 cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 314 370 for q in self.q_vectors 315 371 ] … … 317 373 318 374 def release(self): 375 """ 376 Free the memory. 377 """ 319 378 for b in self.q_buffers: 320 379 b.release() 321 380 self.q_buffers = [] 322 381 382 def __del__(self): 383 self.release() 384 323 385 class GpuKernel(object): 324 386 """ 325 387 Callable SAS kernel. 326 388 327 *kernel* is the GpuKernel object to call .389 *kernel* is the GpuKernel object to call 328 390 329 391 *info* is the module information 330 392 331 *q_input* is the DllInput q vectors at which the kernel should be 332 evaluated. 393 *q_vectors* is the q vectors at which the kernel should be evaluated 394 395 *dtype* is the kernel precision 333 396 334 397 The resulting call method takes the *pars*, a list of values for … … 340 403 Call :meth:`release` when done with the kernel instance. 341 404 """ 342 def __init__(self, kernel, info, q_ input):343 self.q_input = q_input405 def __init__(self, kernel, info, q_vectors, dtype): 406 q_input = GpuInput(q_vectors, dtype) 344 407 self.kernel = kernel 345 408 self.info = info 346 409 self.res = np.empty(q_input.nq, q_input.dtype) 347 dim = '2d' if q_input.is_2 Delse '1d'410 dim = '2d' if q_input.is_2d else '1d' 348 411 self.fixed_pars = info['partype']['fixed-' + dim] 349 412 self.pd_pars = info['partype']['pd-' + dim] … … 352 415 # Note: res may be shorter than res_b if global_size != nq 353 416 env = environment() 354 self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE, 355 2 * MAX_LOOPS * q_input.dtype.itemsize) 356 for _ in env.queues] 357 self.res_b = [cl.Buffer(env.context, mf.READ_WRITE, 358 q_input.global_size[0] * q_input.dtype.itemsize) 359 for _ in env.queues] 360 417 self.queue = env.get_queue(dtype) 418 self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 419 2 * MAX_LOOPS * q_input.dtype.itemsize) 420 self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 421 q_input.global_size[0] * q_input.dtype.itemsize) 422 self.q_input = q_input 423 424 self._need_release = [self.loops_b, self.res_b, self.q_input] 361 425 362 426 def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): … … 366 430 else np.float32) # will never get here, so use np.float32 367 431 368 device_num = 0 369 queuei = environment().queues[device_num] 370 res_bi = self.res_b[device_num] 432 res_bi = self.res_b 371 433 nq = np.uint32(self.q_input.nq) 372 434 if pd_pars: … … 383 445 raise ValueError("too many polydispersity points") 384 446 385 loops_bi = self.loops_b [device_num]386 cl.enqueue_copy( queuei, loops_bi, loops)447 loops_bi = self.loops_b 448 cl.enqueue_copy(self.queue, loops_bi, loops) 387 449 loops_l = cl.LocalMemory(len(loops.data)) 388 450 #ctx = environment().context … … 393 455 fixed = [real(p) for p in fixed_pars] 394 456 args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed 395 self.kernel( queuei, self.q_input.global_size, None, *args)396 cl.enqueue_copy( queuei, self.res, res_bi)457 self.kernel(self.queue, self.q_input.global_size, None, *args) 458 cl.enqueue_copy(self.queue, self.res, res_bi) 397 459 398 460 return self.res 399 461 400 462 def release(self): 401 for b in self.loops_b:402 b.release()403 self.loops_b = []404 for b in self.res_b:405 b.release()406 self. res_b= []463 """ 464 Release resources associated with the kernel. 465 """ 466 for v in self._need_release: 467 v.release() 468 self._need_release = [] 407 469 408 470 def __del__(self): -
sasmodels/kerneldll.py
r74667d3 reafc9fa 1 1 r""" 2 C types wrapper for sasview models. 2 DLL driver for C kernels 3 3 4 4 The global attribute *ALLOW_SINGLE_PRECISION_DLLS* should be set to *True* if … … 44 44 directory for this application, then OpenMP should be supported. 45 45 """ 46 from __future__ import print_function 46 47 47 48 import sys … … 122 123 models are allowed as DLLs. 123 124 """ 124 if callable(info.get('Iq', None)):125 if callable(info.get('Iq', None)): 125 126 return PyModel(info) 126 127 … … 139 140 140 141 source = generate.convert_type(source, dtype) 141 source_files = generate. sources(info) + [info['filename']]142 dll = dll_path(info, dtype)142 source_files = generate.model_sources(info) + [info['filename']] 143 dll = dll_path(info, dtype) 143 144 newest = max(os.path.getmtime(f) for f in source_files) 144 if not os.path.exists(dll) or os.path.getmtime(dll) <newest:145 if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 145 146 # Replace with a proper temp file 146 fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix)147 os.fdopen(fid, "w").write(source)147 fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix) 148 os.fdopen(fid, "w").write(source) 148 149 command = COMPILE%{"source":filename, "output":dll} 149 150 print("Compile command: "+command) … … 160 161 def load_dll(source, info, dtype="double"): 161 162 """ 162 Create and load a dll corresponding to the source, info pair returned163 Create and load a dll corresponding to the source, info pair returned 163 164 from :func:`sasmodels.generate.make` compiled for the target precision. 164 165 … … 199 200 Npd2d = len(self.info['partype']['pd-2d']) 200 201 201 #print("dll", self.dllpath)202 #print("dll", self.dllpath) 202 203 try: 203 204 self.dll = ct.CDLL(self.dllpath) … … 210 211 else c_longdouble) 211 212 pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 212 pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []213 pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 213 214 self.Iq = self.dll[generate.kernel_name(self.info, False)] 214 215 self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d … … 218 219 219 220 def __getstate__(self): 220 return {'info': self.info, 'dllpath': self.dllpath, 'dll': None}221 return self.info, self.dllpath 221 222 222 223 def __setstate__(self, state): 223 self. __dict__= state224 225 def __call__(self, q_input): 226 if self.dtype != q_input.dtype:227 raise TypeError("data is %s kernel is %s" % (q_input.dtype, self.dtype))224 self.info, self.dllpath = state 225 self.dll = None 226 227 def __call__(self, q_vectors): 228 q_input = PyInput(q_vectors, self.dtype) 228 229 if self.dll is None: self._load_dll() 229 kernel = self.Iqxy if q_input.is_2 Delse self.Iq230 kernel = self.Iqxy if q_input.is_2d else self.Iq 230 231 return DllKernel(kernel, self.info, q_input) 231 232 232 # pylint: disable=no-self-use233 def make_input(self, q_vectors):234 """235 Make q input vectors available to the model.236 237 Note that each model needs its own q vector even if the case of238 mixture models because some models may be OpenCL, some may be239 ctypes and some may be pure python.240 """241 return PyInput(q_vectors, dtype=self.dtype)242 243 233 def release(self): 234 """ 235 Release any resources associated with the model. 236 """ 244 237 pass # TODO: should release the dll 245 238 … … 257 250 258 251 The resulting call method takes the *pars*, a list of values for 259 the fixed parameters to the kernel, and *pd_pars*, a list of (value, weight)252 the fixed parameters to the kernel, and *pd_pars*, a list of (value, weight) 260 253 vectors for the polydisperse parameters. *cutoff* determines the 261 254 integration limits: any points with combined weight less than *cutoff* … … 269 262 self.kernel = kernel 270 263 self.res = np.empty(q_input.nq, q_input.dtype) 271 dim = '2d' if q_input.is_2 Delse '1d'264 dim = '2d' if q_input.is_2d else '1d' 272 265 self.fixed_pars = info['partype']['fixed-'+dim] 273 266 self.pd_pars = info['partype']['pd-'+dim] … … 299 292 300 293 def release(self): 294 """ 295 Release any resources associated with the kernel. 296 """ 301 297 pass -
sasmodels/kernelpy.py
r4c2c535 reafc9fa 1 """ 2 Python driver for python kernels 3 4 Calls the kernel with a vector of $q$ values for a single parameter set. 5 Polydispersity is supported by looping over different parameter sets and 6 summing the results. The interface to :class:`PyModel` matches those for 7 :class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`. 8 """ 1 9 import numpy as np 2 10 from numpy import pi, cos … … 5 13 6 14 class PyModel(object): 15 """ 16 Wrapper for pure python models. 17 """ 7 18 def __init__(self, info): 8 19 self.info = info 9 20 10 def __call__(self, q_input): 11 kernel = self.info['Iqxy'] if q_input.is_2D else self.info['Iq'] 21 def __call__(self, q_vectors): 22 q_input = PyInput(q_vectors, dtype=F64) 23 kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq'] 12 24 return PyKernel(kernel, self.info, q_input) 13 25 14 # pylint: disable=no-self-use15 def make_input(self, q_vectors):16 return PyInput(q_vectors, dtype=F64)17 18 26 def release(self): 27 """ 28 Free resources associated with the model. 29 """ 19 30 pass 20 31 … … 41 52 self.nq = q_vectors[0].size 42 53 self.dtype = dtype 43 self.is_2 D= (len(q_vectors) == 2)54 self.is_2d = (len(q_vectors) == 2) 44 55 self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 45 56 self.q_pointers = [q.ctypes.data for q in self.q_vectors] 46 57 47 58 def release(self): 59 """ 60 Free resources associated with the model inputs. 61 """ 48 62 self.q_vectors = [] 49 63 … … 71 85 self.q_input = q_input 72 86 self.res = np.empty(q_input.nq, q_input.dtype) 73 dim = '2d' if q_input.is_2 Delse '1d'87 dim = '2d' if q_input.is_2d else '1d' 74 88 # Loop over q unless user promises that the kernel is vectorized by 75 89 # taggining it with vectorized=True … … 77 91 if dim == '2d': 78 92 def vector_kernel(qx, qy, *args): 93 """ 94 Vectorized 2D kernel. 95 """ 79 96 return np.array([kernel(qxi, qyi, *args) 80 97 for qxi, qyi in zip(qx, qy)]) 81 98 else: 82 99 def vector_kernel(q, *args): 100 """ 101 Vectorized 1D kernel. 102 """ 83 103 return np.array([kernel(qi, *args) 84 104 for qi in q]) … … 122 142 123 143 def release(self): 144 """ 145 Free resources associated with the kernel. 146 """ 124 147 self.q_input = None 125 148 -
sasmodels/model_test.py
r9404dd3 r13ed84c 9 9 if model1 is 'all', then all except the remaining models will be tested 10 10 11 Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1,0.1),11 Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1), 12 12 and the ER and VR are computed. The return values at these points are not 13 13 considered. The test is only to verify that the models run to completion, … … 30 30 [ {parameters}, (qx, qy), I(qx, Iqy)], 31 31 [ {parameters}, [(qx1, qy1), (qx2, qy2), ...], 32 [I(qx1, qy1), I(qx2,qy2), ...]],32 [I(qx1, qy1), I(qx2, qy2), ...]], 33 33 34 34 [ {parameters}, 'ER', ER(pars) ], … … 43 43 Precision defaults to 5 digits (relative). 44 44 """ 45 from __future__ import print_function 45 46 46 47 import sys … … 55 56 56 57 def make_suite(loaders, models): 58 """ 59 Construct the pyunit test suite. 60 61 *loaders* is the list of kernel drivers to use, which is one of 62 *["dll", "opencl"]*, *["dll"]* or *["opencl"]*. For python models, 63 the python driver is always used. 64 65 *models* is the list of models to test, or *["all"]* to test all models. 66 """ 57 67 58 68 ModelTestCase = _hide_model_case_from_nosetests() … … 75 85 # don't try to call cl kernel since it will not be 76 86 # available in some environmentes. 77 is_py = callable(getattr(model_definition, 'Iq', None))87 is_py = callable(getattr(model_definition, 'Iq', None)) 78 88 79 89 if is_py: # kernel implemented in python … … 90 100 test_name = "Model: %s, Kernel: OpenCL"%model_name 91 101 test_method_name = "test_%s_opencl" % model_name 102 # Using dtype=None so that the models that are only 103 # correct for double precision are not tested using 104 # single precision. The choice is determined by the 105 # presence of *single=False* in the model file. 92 106 test = ModelTestCase(test_name, model_definition, 93 107 test_method_name, 94 platform="ocl", dtype= 'single')108 platform="ocl", dtype=None) 95 109 #print("defining", test_name) 96 110 suite.addTest(test) … … 111 125 def _hide_model_case_from_nosetests(): 112 126 class ModelTestCase(unittest.TestCase): 127 """ 128 Test suit for a particular model with a particular kernel driver. 129 130 The test suite runs a simple smoke test to make sure the model 131 functions, then runs the list of tests at the bottom of the model 132 description file. 133 """ 113 134 def __init__(self, test_name, definition, test_method_name, 114 135 platform, dtype): … … 123 144 def _runTest(self): 124 145 smoke_tests = [ 125 [{}, 0.1,None],126 [{}, (0.1,0.1),None],127 [{}, 'ER',None],128 [{}, 'VR',None],146 [{}, 0.1, None], 147 [{}, (0.1, 0.1), None], 148 [{}, 'ER', None], 149 [{}, 'VR', None], 129 150 ] 130 151 … … 140 161 ## values an error. Only do so for the "dll" tests 141 162 ## to reduce noise from both opencl and dll, and because 142 ## python kernels us 163 ## python kernels use platform="dll". 143 164 #raise Exception("No test cases provided") 144 165 pass … … 163 184 actual = [call_VR(model.info, pars)] 164 185 elif isinstance(x[0], tuple): 165 Qx, Qy = zip(*x)186 Qx, Qy = zip(*x) 166 187 q_vectors = [np.array(Qx), np.array(Qy)] 167 188 kernel = make_kernel(model, q_vectors) … … 179 200 # smoke test --- make sure it runs and produces a value 180 201 self.assertTrue(np.isfinite(actual_yi), 181 'invalid f(%s): %s' % (xi, actual_yi))202 'invalid f(%s): %s' % (xi, actual_yi)) 182 203 else: 183 err = abs(yi - actual_yi) 184 nrm = abs(yi) 185 self.assertLess(err * 10**5, nrm, 186 'f(%s); expected:%s; actual:%s' % (xi, yi, actual_yi)) 204 self.assertTrue(is_near(yi, actual_yi, 5), 205 'f(%s); expected:%s; actual:%s' 206 % (xi, yi, actual_yi)) 187 207 188 208 return ModelTestCase 189 209 190 210 def is_near(target, actual, digits=5): 211 """ 212 Returns true if *actual* is within *digits* significant digits of *target*. 213 """ 214 import math 215 shift = 10**math.ceil(math.log10(abs(target))) 216 return abs(target-actual)/shift < 1.5*10**-digits 191 217 192 218 def main(): … … 199 225 200 226 models = sys.argv[1:] 227 if models and models[0] == '-v': 228 verbosity = 2 229 models = models[1:] 230 else: 231 verbosity = 1 201 232 if models and models[0] == 'opencl': 202 233 if not HAVE_OPENCL: … … 217 248 print("""\ 218 249 usage: 219 python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ... 250 python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 251 252 If -v is included on the 253 If neither opencl nor dll is specified, then models will be tested with 254 both opencl and dll; the compute target is ignored for pure python models. 220 255 221 256 If model1 is 'all', then all except the remaining models will be tested. 222 If no compute target is specified, then models will be tested with both opencl 223 and dll; the compute target is ignored for pure python models.""")257 258 """) 224 259 225 260 return 1 226 261 227 262 #runner = unittest.TextTestRunner() 228 runner = xmlrunner.XMLTestRunner(output='logs' )263 runner = xmlrunner.XMLTestRunner(output='logs', verbosity=verbosity) 229 264 result = runner.run(make_suite(loaders, models)) 230 265 return 1 if result.failures or result.errors else 0 … … 237 272 Run "nosetests sasmodels" on the command line to invoke it. 238 273 """ 239 tests = make_suite(['opencl', 'dll'],['all'])274 tests = make_suite(['opencl', 'dll'], ['all']) 240 275 for test_i in tests: 241 276 yield test_i._runTest -
sasmodels/models/HayterMSAsq.py
reb69cce r13ed84c 56 56 parameters used in P(Q). 57 57 """ 58 single = False # double precision only for now 58 59 # [ "name", "units", default, [lower, upper], "type", "description" ], 59 60 parameters = [["effect_radius", "Ang", 20.75, [0, inf], "volume", … … 90 91 # odd that the default st has saltconc zero 91 92 demo = dict(effect_radius = 20.75,charge=19.0,volfraction = 0.0192,temperature=318.16,saltconc=0.05,dielectconst=71.08,effect_radius_pd = 0.1,effect_radius_pd_n = 40) 93 # 94 # attempt to use same values as old sasview unit test 95 tests = [ 96 [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 20.75, 'charge' : 19.0, 'volfraction' : 0.0192, 'temperature' : 298.0, 97 'saltconc' : 0,'dielectconst' : 78.0, 'effect_radius_pd' : 0}, [0.0010], [0.0712928]] 98 ] 92 99 -
sasmodels/models/barbell.py
reb69cce rdcdf29d 100 100 title = "Cylinder with spherical end caps" 101 101 description = """ 102 Calculates the scattering from a barbell-shaped cylinder. That is a sphereocylinder with spherical end caps that have a radius larger than that of the cylinder and the center of the end cap 103 radius lies outside of the cylinder. 104 Note: As the length of cylinder(bar) -->0,it becomes a dumbbell. And when rad_bar = rad_bell, it is a spherocylinder. 105 It must be that rad_bar <(=) rad_bell. 102 Calculates the scattering from a barbell-shaped cylinder. 103 That is a sphereocylinder with spherical end caps that have a radius larger 104 than that of the cylinder and the center of the end cap radius lies outside 105 of the cylinder. 106 Note: As the length of cylinder(bar) -->0,it becomes a dumbbell. And when 107 rad_bar = rad_bell, it is a spherocylinder. 108 It must be that rad_bar <(=) rad_bell. 106 109 """ 107 110 category = "shape:cylinder" 108 111 # pylint: disable=bad-whitespace, line-too-long 109 112 # ["name", "units", default, [lower, upper], "type","description"], 110 parameters = [["sld", "4e-6/Ang^2", 4, [-inf, inf], "","Barbell scattering length density"],111 ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "","Solvent scattering length density"],112 ["bell_radius", "Ang", 40, [0, inf], "volume","Spherical bell radius"],113 ["radius", "Ang", 20, [0, inf], "volume","Cylindrical bar radius"],114 ["length", "Ang", 400, [0, inf], "volume","Cylinder bar length"],115 ["theta", "degrees",60, [-inf, inf], "orientation", "In plane angle"],116 ["phi", "degrees",60, [-inf, inf], "orientation", "Out of plane angle"],113 parameters = [["sld", "4e-6/Ang^2", 4, [-inf, inf], "", "Barbell scattering length density"], 114 ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"], 115 ["bell_radius", "Ang", 40, [0, inf], "volume", "Spherical bell radius"], 116 ["radius", "Ang", 20, [0, inf], "volume", "Cylindrical bar radius"], 117 ["length", "Ang", 400, [0, inf], "volume", "Cylinder bar length"], 118 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 119 ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 117 120 ] 121 # pylint: enable=bad-whitespace, line-too-long 118 122 119 123 source = ["lib/J1.c", "lib/gauss76.c", "barbell.c"] -
sasmodels/models/bcc.py
reb69cce r13ed84c 109 109 title = "Body-centred cubic lattic with paracrystalline distortion" 110 110 description = """ 111 Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations 112 are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is 113 assumed to be isotropic and characterized by a Gaussian distribution. 111 Calculates the scattering from a **body-centered cubic lattice** with 112 paracrystalline distortion. Thermal vibrations are considered to be 113 negligible, and the size of the paracrystal is infinitely large. 114 Paracrystalline distortion is assumed to be isotropic and characterized 115 by a Gaussian distribution. 114 116 """ 115 117 category = "shape:paracrystal" 116 118 119 single = False 120 121 # pylint: disable=bad-whitespace, line-too-long 117 122 # ["name", "units", default, [lower, upper], "type","description" ], 118 parameters = [["dnn", "Ang", 220, [-inf, inf], "","Nearest neighbour distance"],119 ["d_factor", "", 0.06, [-inf, inf], "","Paracrystal distortion factor"],120 ["radius", "Ang", 40, [0, inf], "volume","Particle radius"],121 ["sld", "1e-6/Ang^2", 4, [-inf, inf], "","Particle scattering length density"],122 ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "","Solvent scattering length density"],123 ["theta", "degrees", 60,[-inf, inf], "orientation", "In plane angle"],124 ["phi", "degrees", 60,[-inf, inf], "orientation", "Out of plane angle"],125 ["psi", "degrees", 60,[-inf, inf], "orientation", "Out of plane angle"]123 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"], 124 ["d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"], 125 ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"], 126 ["sld", "1e-6/Ang^2", 4, [-inf, inf], "", "Particle scattering length density"], 127 ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"], 128 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 129 ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 130 ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"] 126 131 ] 132 # pylint: enable=bad-whitespace, line-too-long 127 133 128 134 source = ["lib/J1.c", "lib/gauss150.c", "bcc.c"] -
sasmodels/models/broad_peak.py
reb69cce rdcdf29d 60 60 category = "shape-independent" 61 61 62 # pylint: disable=bad-whitespace, line-too-long 62 63 # ["name", "units", default, [lower, upper], "type", "description"], 63 parameters = [["porod_scale", "",1.0e-05, [-inf, inf], "", "Power law scale factor"],64 ["porod_exp", "",3.0, [-inf, inf], "", "Exponent of power law"],65 ["lorentz_scale", "",10.0, [-inf, inf], "", "Scale factor for broad Lorentzian peak"],66 ["lorentz_length", "Ang", 50.0, [-inf, inf], "", "Lorentzian screening length"],67 ["peak_pos", "1/Ang", 0.1, [-inf, inf], "", "Peak postion in q"],68 ["lorentz_exp", "", 2.0, [-inf, inf], "", "exponent of Lorentz function"],64 parameters = [["porod_scale", "", 1.0e-05, [-inf, inf], "", "Power law scale factor"], 65 ["porod_exp", "", 3.0, [-inf, inf], "", "Exponent of power law"], 66 ["lorentz_scale", "", 10.0, [-inf, inf], "", "Scale factor for broad Lorentzian peak"], 67 ["lorentz_length", "Ang", 50.0, [-inf, inf], "", "Lorentzian screening length"], 68 ["peak_pos", "1/Ang", 0.1, [-inf, inf], "", "Peak position in q"], 69 ["lorentz_exp", "", 2.0, [-inf, inf], "", "Exponent of Lorentz function"], 69 70 ] 71 # pylint: enable=bad-whitespace, line-too-long 70 72 71 def Iq(q, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp): 73 def Iq(q, 74 porod_scale=1.0e-5, 75 porod_exp=3.0, 76 lorentz_scale=10.0, 77 lorentz_length=50.0, 78 peak_pos=0.1, 79 lorentz_exp=2.0): 80 """ 81 :param q: Input q-value 82 :param porod_scale: Power law scale factor 83 :param porod_exp: Exponent of power law 84 :param lorentz_scale: Scale factor for broad Lorentzian peak 85 :param lorentz_length: Lorentzian screening length 86 :param peak_pos: Peak position in q 87 :param lorentz_exp: Exponent of Lorentz function 88 :return: Calculated intensity 89 """ 90 72 91 inten = (porod_scale / q ** porod_exp + lorentz_scale 73 92 / (1.0 + (abs(q - peak_pos) * lorentz_length) ** lorentz_exp)) … … 76 95 77 96 def Iqxy(qx, qy, *args): 97 """ 98 :param qx: Input q_x-value 99 :param qy: Input q_y-value 100 :param args: Remaining arguments 101 :return: 2D-Intensity 102 """ 78 103 return Iq(sqrt(qx ** 2 + qy ** 2), *args) 104 79 105 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 80 81 106 82 107 demo = dict(scale=1, background=0, -
sasmodels/models/capped_cylinder.py
reb69cce rdcdf29d 125 125 """ 126 126 category = "shape:cylinder" 127 # pylint: disable=bad-whitespace, line-too-long 128 # ["name", "units", default, [lower, upper], "type", "description"], 129 parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "", "Cylinder scattering length density"], 130 ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", "Solvent scattering length density"], 131 ["radius", "Ang", 20, [0, inf], "volume", "Cylinder radius"], 127 132 128 # ["name", "units", default, [lower, upper], "type", "description"],129 parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "",130 "Cylinder scattering length density"],131 ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",132 "Solvent scattering length density"],133 ["radius", "Ang", 20, [0, inf], "volume", "Cylinder radius"],134 133 # TODO: use an expression for cap radius with fixed bounds. 135 134 # The current form requires cap radius R bigger than cylinder radius r. … … 141 140 # in the capped cylinder, and zero for no bar in the barbell model. In 142 141 # both models, one would be a pill. 143 ["cap_radius", "Ang", 20, [0, inf],"volume", "Cap radius"],144 ["length", "Ang", 400, [0, inf],"volume", "Cylinder length"],145 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"],146 ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"],142 ["cap_radius", "Ang", 20, [0, inf], "volume", "Cap radius"], 143 ["length", "Ang", 400, [0, inf], "volume", "Cylinder length"], 144 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 145 ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 147 146 ] 147 # pylint: enable=bad-whitespace, line-too-long 148 148 149 149 source = ["lib/J1.c", "lib/gauss76.c", "capped_cylinder.c"] … … 159 159 phi_pd=15, phi_pd_n=1) 160 160 oldname = 'CappedCylinderModel' 161 oldpars = dict(sld='sld_capcyl', solvent_sld='sld_solv', 162 length='len_cyl', radius='rad_cyl', cap_radius='rad_cap') 161 oldpars = dict(sld='sld_capcyl', 162 solvent_sld='sld_solv', 163 length='len_cyl', 164 radius='rad_cyl', 165 cap_radius='rad_cap') -
sasmodels/models/fcc.c
r82d239a reeb8bac 101 101 102 102 double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 103 double q_z;103 // double q_z; 104 104 double cos_val_b3, cos_val_b2, cos_val_b1; 105 105 double a1_dot_q, a2_dot_q,a3_dot_q; … … 124 124 const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1); 125 125 // q vector 126 q_z = 0.0; // for SANS; assuming qz is negligible126 // q_z = 0.0; // for SANS; assuming qz is negligible 127 127 /// Angles here are respect to detector coordinate 128 128 /// instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) -
sasmodels/models/fcc.py
reb69cce r13ed84c 112 112 category = "shape:paracrystal" 113 113 114 single = False 115 114 116 # ["name", "units", default, [lower, upper], "type","description"], 115 117 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"], -
sasmodels/models/gaussian_peak.py
reb69cce r13ed84c 42 42 category = "shape-independent" 43 43 44 single = False 44 45 # ["name", "units", default, [lower, upper], "type","description"], 45 46 parameters = [["q0", "1/Ang", 0.05, [-inf, inf], "", "Peak position"], -
sasmodels/models/gel_fit.py
r513efc5 r168052c 7 7 in the position of the polymer chains that ensure thermodynamic equilibrium, 8 8 and a longer distance (denoted here as $a2$ ) needed to account for the static 9 accumulations of polymer pinned down by junction points or clusters of such points.10 The latter is derived from a simple Guinier function.9 accumulations of polymer pinned down by junction points or clusters of such 10 points. The latter is derived from a simple Guinier function. 11 11 12 12 … … 40 40 --------- 41 41 42 Mitsuhiro Shibayama, Toyoichi Tanaka, Charles C Han, J. Chem. Phys. 1992, 97 (9),43 6829-684142 Mitsuhiro Shibayama, Toyoichi Tanaka, Charles C Han, 43 *J. Chem. Phys.* 1992, 97 (9), 6829-6841 44 44 45 45 Simon Mallam, Ferenc Horkay, Anne-Marie Hecht, Adrian R Rennie, Erik Geissler, 46 Macromolecules1991, 24, 543-54846 *Macromolecules* 1991, 24, 543-548 47 47 48 48 """ … … 62 62 category = "shape-independent" 63 63 64 # pylint: disable=bad-whitespace, line-too-long 64 65 # ["name", "units", default, [lower, upper], "type","description"], 65 66 parameters = [["guinier_scale", "cm^{-1}", 1.7, [-inf, inf], "", "Guinier length scale"], … … 68 69 ["fractal_exp", "", 2.0, [0, inf], "", "Fractal exponent"], 69 70 ["cor_length", "Ang", 16.0, [0, inf], "", "Correlation length"] 70 71 71 ] 72 # pylint: enable=bad-whitespace, line-too-long 72 73 73 74 source = ["gel_fit.c"] … … 92 93 'fractal_exp': 10.0, 93 94 'cor_length': 20.0 94 95 }, 0.1, 0.716532], 95 96 96 97 [{'guinier_scale': 4.0, … … 100 101 'cor_length': 20.0, 101 102 'background': 20.0, 102 }, 5.0, 20.1224653026], 103 104 ] 103 }, 5.0, 20.1224653026], 104 ] -
sasmodels/models/hardsphere.py
r7ed702f r13ed84c 55 55 "volume fraction of hard spheres"], 56 56 ] 57 single = False 57 58 58 59 # No volume normalization despite having a volume parameter … … 106 107 oldpars = dict() 107 108 109 tests = [ 110 [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 50.0, 'volfraction' : 0.2, 111 'effect_radius_pd' : 0}, [0.001], [0.209128]] 112 ] 113 -
sasmodels/models/lamellarCaille.py
rd18f8a8 r13ed84c 87 87 category = "shape:lamellae" 88 88 89 single = False 90 89 91 # ["name", "units", default, [lower, upper], "type","description"], 90 92 parameters = [["thickness", "Ang", 30.0, [0, inf], "volume", "sheet thickness"], … … 122 124 oldpars = dict(thickness='delta', Nlayers='N_plates', Caille_parameter='caille', 123 125 sld='sld_bi',solvent_sld='sld_sol') 126 # 127 tests = [ 128 [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 30.,'Nlayers' : 20.0, 'spacing' : 400., 129 'Caille_parameter' : 0.1, 'sld' : 6.3, 'solvent_sld' : 1.0, 130 'thickness_pd' : 0.0, 'spacing_pd' : 0.0 }, [0.001], [28895.13397]] 131 ] -
sasmodels/models/lamellarCailleHG.py
rd18f8a8 r13ed84c 91 91 category = "shape:lamellae" 92 92 93 single = False 93 94 parameters = [ 94 95 96 [ "tail_length", "Ang",10, [0, inf], "volume",97 "Tail thickness"],98 [ "head_length", "Ang",2, [0, inf], "volume",99 "head thickness"],100 [ "Nlayers", "",30, [0, inf], "",101 "Number of layers"],102 [ "spacing", "Ang", 40., [0.0,inf], "volume",103 "d-spacing of Caille S(Q)"],104 [ "Caille_parameter", "", 0.001, [0.0,0.8], "",105 "Caille parameter"],106 [ "sld", "1e-6/Ang^2", 0.4, [-inf,inf], "",107 "Tail scattering length density"],108 [ "head_sld", "1e-6/Ang^2", 2.0, [-inf,inf], "",109 "Head scattering length density"],110 [ "solvent_sld", "1e-6/Ang^2", 6, [-inf,inf], "",111 "Solvent scattering length density"],95 # [ "name", "units", default, [lower, upper], "type", 96 # "description" ], 97 ["tail_length", "Ang", 10, [0, inf], "volume", 98 "Tail thickness"], 99 ["head_length", "Ang", 2, [0, inf], "volume", 100 "head thickness"], 101 ["Nlayers", "", 30, [0, inf], "", 102 "Number of layers"], 103 ["spacing", "Ang", 40., [0.0, inf], "volume", 104 "d-spacing of Caille S(Q)"], 105 ["Caille_parameter", "", 0.001, [0.0, 0.8], "", 106 "Caille parameter"], 107 ["sld", "1e-6/Ang^2", 0.4, [-inf, inf], "", 108 "Tail scattering length density"], 109 ["head_sld", "1e-6/Ang^2", 2.0, [-inf, inf], "", 110 "Head scattering length density"], 111 ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "", 112 "Solvent scattering length density"], 112 113 ] 113 114 114 source = [ 115 source = ["lamellarCailleHG_kernel.c"] 115 116 116 117 # No volume normalization despite having a volume parameter … … 128 129 129 130 demo = dict( 130 scale=1, background=0, 131 Nlayers=20, 132 spacing=200., Caille_parameter=0.05, 133 tail_length=15,head_length=10, 134 #sld=-1, head_sld=4.0, solvent_sld=6.0, 135 sld=-1, head_sld=4.1, solvent_sld=6.0, 136 tail_length_pd= 0.1, tail_length_pd_n=20, 137 head_length_pd= 0.05, head_length_pd_n=30, 138 spacing_pd= 0.2, spacing_pd_n=40 139 ) 131 scale=1, background=0, 132 Nlayers=20, spacing=200., Caille_parameter=0.05, 133 tail_length=15, head_length=10, 134 #sld=-1, head_sld=4.0, solvent_sld=6.0, 135 sld=-1, head_sld=4.1, solvent_sld=6.0, 136 tail_length_pd=0.1, tail_length_pd_n=20, 137 head_length_pd=0.05, head_length_pd_n=30, 138 spacing_pd=0.2, spacing_pd_n=40, 139 ) 140 140 141 141 oldname = 'LamellarPSHGModel' 142 oldpars = dict(tail_length='deltaT',head_length='deltaH',Nlayers='n_plates',Caille_parameter='caille', sld='sld_tail', head_sld='sld_head',solvent_sld='sld_solvent') 142 143 oldpars = dict( 144 tail_length='deltaT', head_length='deltaH', Nlayers='n_plates', 145 Caille_parameter='caille', sld='sld_tail', head_sld='sld_head', 146 solvent_sld='sld_solvent') 147 # 148 tests = [ 149 [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 10.0, 'head_length' : 2.0,'Nlayers' : 30.0, 'spacing' : 40., 150 'Caille_parameter' : 0.001, 'sld' : 0.4, 'head_sld' : 2.0, 'solvent_sld' : 6.0, 151 'tail_length_pd' : 0.0, 'head_length_pd' : 0.0, 'spacing_pd' : 0.0 }, [0.001], [6838238.571488]] 152 ] -
sasmodels/models/lamellarFFHG.py
reb69cce r7f47777 121 121 oldpars = dict(head_length='h_thickness', tail_length='t_length', 122 122 sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent') 123 # 124 tests = [ 125 [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 15.0, 'head_length' : 10.0,'sld' : 0.4, 126 'head_sld' : 3.0, 'solvent_sld' : 6.0, }, [0.001], [653143.9209]] 127 ] 123 128 129 -
sasmodels/models/lamellarPC.py
rd18f8a8 r13ed84c 111 111 category = "shape:lamellae" 112 112 113 single = False 114 113 115 # ["name", "units", default, [lower, upper], "type","description"], 114 116 parameters = [["thickness", "Ang", 33.0, [0, inf], "volume", … … 148 150 oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer', 149 151 solvent_sld='sld_solvent') 152 # 153 tests = [ 154 [ {'scale': 1.0, 'background' : 0.0, 'thickness' : 33.,'Nlayers' : 20.0, 'spacing' : 250., 'spacing_polydisp' : 0.2, 155 'sld' : 1.0, 'solvent_sld' : 6.34, 156 'thickness_pd' : 0.0, 'thickness_pd_n' : 40 }, [0.001, 0.215268], [21829.3, 0.00487686]] 157 ] -
sasmodels/models/lib/Si.c
r7d256c8 rf12357f 13 13 14 14 // Explicitly writing factorial values triples the speed of the calculation 15 out_cos = 1 /x - 2/pow(x,3) + 24/pow(x,5) - 720/pow(x,7) + 40320/pow(x,9);16 out_sin = 1 /x - 6/pow(x,4) + 120/pow(x,6) - 5040/pow(x,8) + 362880/pow(x,10);15 out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7); 16 out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8); 17 17 18 18 out -= cos(x) * out_cos; … … 22 22 23 23 // Explicitly writing factorial values triples the speed of the calculation 24 out = x - pow(x, 3)/18 + pow(x,5)/600 - pow(x,7)/35280 + pow(x,9)/3265920;24 out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.; 25 25 26 //printf ("Si=%g %g\n", x, out);27 26 return out; 28 27 } -
sasmodels/models/lib/sph_j1c.c
r9c461c7 ra3e78c3 1 1 /** 2 * Spherical Bessel function j1(x)/x2 * Spherical Bessel function 3*j1(x)/x 3 3 * 4 4 * Used for low q to avoid cancellation error. … … 8 8 * requires the first 3 terms. Double precision requires the 4th term. 9 9 * The fifth term is not needed, and is commented out. 10 * Taylor expansion: 11 * 1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.)))) 12 * Expression returned from Herbie (herbie.uwpise.org/demo): 13 * const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); 14 * return t*t; 10 15 */ 11 16 … … 18 23 SINCOS(q, sin_q, cos_q); 19 24 20 const double bessel = (q < 1.e-1) ?21 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))) // + q2*(3./3991680.)))25 const double bessel = (q < 0.384038453352533) 26 ? (1.0 + q2*(-3./30. + q2*(3./840.))) 22 27 : 3.0*(sin_q/q - cos_q)/q2; 23 28 24 29 return bessel; 30 31 /* 32 // Code to test various expressions 33 if (sizeof(q2) > 4) { 34 return 3.0*(sin_q/q - cos_q)/q2; 35 } else if (q < 0.384038453352533) { 36 //const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); return t*t; 37 return 1.0 + q2*q2*(3./840.) - q2*(3./30.); 38 //return 1.0 + q2*(-3./30. + q2*(3./840.)); 39 //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); 40 //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); 41 } else { 42 return 3.0*(sin_q/q - cos_q)/q2; 43 } 44 */ 25 45 } -
sasmodels/models/mass_fractal.py
r9c461c7 r168052c 52 52 --------- 53 53 54 D Mildner and P Hall, *J. Phys. D: Appl. Phys.*, 19 (1986) 1535-1545 Equation(9) 54 D Mildner and P Hall, *J. Phys. D: Appl. Phys.*, 55 19 (1986) 1535-1545 Equation(9) 55 56 56 57 … … 79 80 category = "shape-independent" 80 81 82 # pylint: disable=bad-whitespace, line-too-long 81 83 # ["name", "units", default, [lower, upper], "type","description"], 82 84 parameters = [["radius", "Ang", 10.0, [0.0, inf], "", "Particle radius"], 83 85 ["mass_dim", "", 1.9, [1.0, 6.0], "", "Mass fractal dimension"], 84 86 ["cutoff_length", "Ang", 100.0, [0.0, inf], "", "Cut-off length"], 85 86 87 ] 88 # pylint: enable=bad-whitespace, line-too-long 87 89 88 90 source = ["lib/sph_j1c.c", "lib/lanczos_gamma.c", "mass_fractal.c"] … … 91 93 radius=10.0, 92 94 mass_dim=1.9, 93 cutoff_length= 2.3)95 cutoff_length=100.0) 94 96 95 97 oldname = 'MassFractalModel' … … 98 100 cutoff_length='co_length') 99 101 100 tests = [[{'radius': 2.0, 101 'mass_dim': 3.3, 102 'cutoff_length': 1.0, 103 }, 0.5, 1.29016774904], 102 tests = [ 104 103 105 [{'radius': 1.0,106 'mass_dim': 1.3,107 'cutoff_length': 1.0,108 'background': 0.8,109 }, 0.001, 1.69747015932],104 # Accuracy tests based on content in test/utest_other_models.py 105 [{'radius': 10.0, 106 'mass_dim': 1.9, 107 'cutoff_length': 100.0, 108 }, 0.05, 279.59322], 110 109 111 [{'radius': 1.0, 112 'mass_dim': 2.3, 113 'cutoff_length': 1.0, 114 'scale': 10.0, 115 }, 0.051, 11.6227966145], 116 ] 110 # Additional tests with larger range of parameters 111 [{'radius': 2.0, 112 'mass_dim': 3.3, 113 'cutoff_length': 1.0, 114 }, 0.5, 1.29016774904], 115 116 [{'radius': 1.0, 117 'mass_dim': 1.3, 118 'cutoff_length': 1.0, 119 'background': 0.8, 120 }, 0.001, 1.69747015932], 121 122 [{'radius': 1.0, 123 'mass_dim': 2.3, 124 'cutoff_length': 1.0, 125 'scale': 10.0, 126 }, 0.051, 11.6227966145], 127 ] -
sasmodels/models/mass_surface_fractal.py
r9c461c7 r168052c 54 54 P Schmidt, *J Appl. Cryst.*, 24 (1991) 414-435 Equation(19) 55 55 56 A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 35 (1987) 2361-2364 Equation(2) 56 A J Hurd, D W Schaefer, J E Martin, *Phys. Rev. A*, 57 35 (1987) 2361-2364 Equation(2) 57 58 58 59 """ … … 79 80 category = "shape-independent" 80 81 82 # pylint: disable=bad-whitespace, line-too-long 81 83 # ["name", "units", default, [lower, upper], "type","description"], 82 84 parameters = [["mass_dim", "", 1.8, [1e-16, 6.0], "", … … 88 90 ["primary_rg", "Ang", 4000., [0.0, inf], "", 89 91 "Primary particle radius of gyration"], 90 91 92 ] 93 # pylint: enable=bad-whitespace, line-too-long 92 94 93 95 source = ["mass_surface_fractal.c"] … … 106 108 primary_rg='primary_rg') 107 109 108 tests = [[{'radius': 2.0, 109 'mass_dim': 3.3, 110 'surface_dim': 1.0, 111 'cluster_rg': 90.0, 112 'primary_rg': 4000.0, 113 }, 0.001, 0.18462699016], 110 tests = [ 114 111 115 [{'radius': 1.0, 116 'mass_dim': 1.3, 117 'surface_dim': 1.0, 118 'cluster_rg': 90.0, 119 'primary_rg': 2000.0, 120 'background': 0.8, 121 }, 0.001, 1.16539753641], 112 # Accuracy tests based on content in test/utest_other_models.py 113 [{'mass_dim': 1.8, 114 'surface_dim': 2.3, 115 'cluster_rg': 86.7, 116 'primary_rg': 4000.0, 117 }, 0.05, 1.77537e-05], 122 118 123 [{'radius': 1.0, 124 'mass_dim': 2.3, 125 'surface_dim': 1.0, 126 'cluster_rg': 90.0, 127 'primary_rg': 1000.0, 128 'scale': 10.0, 129 }, 0.051, 0.000169548800377], 130 ] 119 # Additional tests with larger range of parameters 120 [{'mass_dim': 3.3, 121 'surface_dim': 1.0, 122 'cluster_rg': 90.0, 123 'primary_rg': 4000.0, 124 }, 0.001, 0.18462699016], 125 126 [{'mass_dim': 1.3, 127 'surface_dim': 1.0, 128 'cluster_rg': 90.0, 129 'primary_rg': 2000.0, 130 'background': 0.8, 131 }, 0.001, 1.16539753641], 132 133 [{'mass_dim': 2.3, 134 'surface_dim': 1.0, 135 'cluster_rg': 90.0, 136 'primary_rg': 1000.0, 137 'scale': 10.0, 138 }, 0.051, 0.000169548800377], 139 ] -
sasmodels/models/pearl_necklace.c
rcf404cb rf12357f 12 12 double string_thickness, double number_of_pearls, double sld, 13 13 double string_sld, double solvent_sld); 14 15 double ER(double radius, double edge_separation,16 double string_thickness, double number_of_pearls);17 double VR(double radius, double edge_separation,18 double string_thickness, double number_of_pearls);19 14 20 15 // From Igor library … … 22 17 double num_pearls, double sld_pearl, double sld_string, double sld_solv) 23 18 { 19 //relative slds 24 20 double contrast_pearl = sld_pearl - sld_solv; 25 21 double contrast_string = sld_string - sld_solv; 22 23 // number of string segments 24 num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls 25 double num_strings = num_pearls - 1.0; 26 27 //Pi 28 double pi = 4.0*atan(1.0); 29 30 // center to center distance between the neighboring pearls 31 double A_s = edge_separation + 2.0 * radius; 32 33 // Repeated Calculations 34 double sincasq = sinc(q*A_s); 35 double oneminussinc = 1 - sincasq; 36 double q_r = q * radius; 37 double q_edge = q * edge_separation; 38 39 // each volume 40 double string_vol = edge_separation * pi * thick_string * thick_string / 4.0; 41 double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 26 42 27 43 //total volume 28 double pi = 4.0*atan(1.0); 29 double tot_vol = form_volume(radius, edge_separation, thick_string, num_pearls); 30 double string_vol = edge_separation * pi * pow((thick_string / 2.0), 2); 31 double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3); 32 double num_strings = num_pearls - 1; 33 //each mass 44 double tot_vol; 45 //each masses 34 46 double m_r= contrast_string * string_vol; 35 47 double m_s= contrast_pearl * pearl_vol; 36 48 double psi, gamma, beta; 37 49 //form factors 38 double sss; //pearls 39 double srr; //strings 40 double srs; //cross 41 double A_s, srr_1, srr_2, srr_3; 50 double sss, srr, srs; //cross 51 double srr_1, srr_2, srr_3; 42 52 double form_factor; 53 tot_vol = num_strings * string_vol; 54 tot_vol += num_pearls * pearl_vol; 43 55 44 56 //sine functions of a pearl 45 psi = sin(q*radius); 46 psi -= q * radius * cos(q * radius); 47 psi /= pow((q * radius), 3); 57 psi = sin(q_r); 58 psi -= q_r * cos(q_r); 48 59 psi *= 3.0; 60 psi /= q_r * q_r * q_r; 49 61 50 62 // Note take only 20 terms in Si series: 10 terms may be enough though. 51 gamma = Si(q* edge_separation); 52 gamma /= (q* edge_separation); 53 beta = Si(q * (edge_separation + radius)); 54 beta -= Si(q * radius); 55 beta /= (q* edge_separation); 56 57 // center to center distance between the neighboring pearls 58 A_s = edge_separation + 2.0 * radius; 63 gamma = Si(q_edge); 64 gamma /= (q_edge); 65 beta = Si(q * (A_s - radius)); 66 beta -= Si(q_r); 67 beta /= q_edge; 59 68 60 69 // form factor for num_pearls 61 sss = 1.0 - pow(sinc (q*A_s), num_pearls);62 sss /= pow((1.0-sinc(q*A_s)), 2);63 sss *= -sinc (q*A_s);64 sss -= num_pearls /2.0;65 sss += num_pearls /(1.0-sinc(q*A_s));66 sss *= 2.0 * pow((m_s*psi), 2);70 sss = 1.0 - pow(sincasq, num_pearls); 71 sss /= oneminussinc * oneminussinc; 72 sss *= -sincasq; 73 sss -= num_pearls / 2.0; 74 sss += num_pearls / oneminussinc; 75 sss *= 2.0 * m_s * psi * m_s * psi; 67 76 68 77 // form factor for num_strings (like thin rods) 69 srr_1 = - pow(sinc(q*edge_separation/2.0), 2);78 srr_1 = -sinc(q_edge/2.0) * sinc(q_edge/2.0); 70 79 71 80 srr_1 += 2.0 * gamma; 72 81 srr_1 *= num_strings; 73 srr_2 = 2.0/(1.0-sinc(q*A_s)); 74 srr_2 *= num_strings * pow(beta, 2); 75 srr_3 = 1.0 - pow(sinc(q*A_s), num_strings); 76 srr_3 /= pow((1.0-sinc(q*A_s)), 2); 77 srr_3 *= -2.0 * pow(beta, 2); 82 srr_2 = 2.0/oneminussinc; 83 srr_2 *= num_strings; 84 srr_2 *= beta * beta; 85 srr_3 = 1.0 - pow(sincasq, num_strings); 86 srr_3 /= oneminussinc * oneminussinc; 87 srr_3 *= beta * beta; 88 srr_3 *= -2.0; 78 89 79 90 // total srr 80 91 srr = srr_1 + srr_2 + srr_3; 81 srr *= pow(m_r, 2);92 srr *= m_r * m_r; 82 93 83 94 // form factor for correlations 84 95 srs = 1.0; 85 srs -= pow(sinc(q*A_s), num_strings); 86 srs /= pow((1.0-sinc(q*A_s)), 2); 87 srs *= -sinc(q*A_s); 88 srs += (num_strings/(1.0-sinc(q*A_s))); 89 srs *= 4.0 * (m_r * m_s * beta * psi); 96 srs -= pow(sincasq, num_strings); 97 srs /= oneminussinc * oneminussinc; 98 srs *= -sincasq; 99 srs += num_strings/oneminussinc; 100 srs *= 4.0; 101 srs *= (m_r * m_s * beta * psi); 90 102 91 103 form_factor = sss + srr + srs; … … 103 115 double number_of_strings = number_of_pearls - 1.0; 104 116 105 double string_vol = edge_separation * pi * pow((string_thickness / 2.0), 2);106 double pearl_vol = 4.0 / 3.0 * pi * pow(radius, 3);117 double string_vol = edge_separation * pi * string_thickness * string_thickness / 4.0; 118 double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius; 107 119 108 120 total_vol = number_of_strings * string_vol; … … 124 136 double string_sld, double solvent_sld) 125 137 { 126 double value = 0.0; 127 double tot_vol = 0.0; 138 double value, tot_vol; 139 140 if (string_thickness >= radius || number_of_pearls <= 0) { 141 return NAN; 142 } 128 143 129 144 value = _pearl_necklace_kernel(q, radius, edge_separation, string_thickness, … … 142 157 sld, string_sld, solvent_sld)); 143 158 } 144 145 146 double ER(double radius, double edge_separation,147 double string_thickness, double number_of_pearls)148 {149 double tot_vol = form_volume(radius, edge_separation, string_thickness, number_of_pearls);150 double pi = 4.0*atan(1.0);151 152 double rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333);153 154 return rad_out;155 }156 double VR(double radius, double edge_separation,157 double string_thickness, double number_of_pearls)158 {159 return 1.0;160 } -
sasmodels/models/pearl_necklace.py
rcf404cb rd18582e 1 1 r""" 2 This model provides the form factor for a pearl necklace composed of two 3 elements: *N* pearls (homogeneous spheres of radius *R*) freely jointed by *M* 2 This model provides the form factor for a pearl necklace composed of two 3 elements: *N* pearls (homogeneous spheres of radius *R*) freely jointed by *M* 4 4 rods (like strings - with a total mass *Mw* = *M* \* *m*\ :sub:`r` + *N* \* *m*\ 5 :sub:`s`, and the string segment length (or edge separation) *l* 5 :sub:`s`, and the string segment length (or edge separation) *l* 6 6 (= *A* - 2\ *R*)). *A* is the center-to-center pearl separation distance. 7 7 … … 13 13 ---------- 14 14 15 The output of the scattering intensity function for the PearlNecklaceModel is 15 The output of the scattering intensity function for the PearlNecklaceModel is 16 16 given by (Schweins, 2004) 17 17 … … 37 37 \beta(q) &= \frac{\int_{qR}^{q(A-R)}\frac{sin(t)}{t}dt}{ql} 38 38 39 where the mass *m*\ :sub:`i` is (SLD\ :sub:`i` - SLD\ :sub:`solvent`) \* 39 where the mass *m*\ :sub:`i` is (SLD\ :sub:`i` - SLD\ :sub:`solvent`) \* 40 40 (volume of the *N* pearls/rods). *V* is the total volume of the necklace. 41 41 42 The 2D scattering intensity is the same as $P(q)$ above, regardless of the 42 The 2D scattering intensity is the same as $P(q)$ above, regardless of the 43 43 orientation of the *q* vector. 44 44 … … 54 54 REFERENCE 55 55 56 R Schweins and K Huber, *Particle Scattering Factor of Pearl Necklace Chains*, 56 R Schweins and K Huber, *Particle Scattering Factor of Pearl Necklace Chains*, 57 57 *Macromol. Symp.* 211 (2004) 25-42 2004 58 58 """ 59 59 60 from numpy import inf 60 from numpy import inf, pi 61 61 62 62 name = "pearl_necklace" … … 79 79 80 80 # ["name", "units", default, [lower, upper], "type","description"], 81 parameters = [["radius", "Angstrom", 80.0, [0, inf], "volume", 81 parameters = [["radius", "Angstrom", 80.0, [0, inf], "volume", 82 82 "Mean radius of the chained spheres"], 83 ["edge_separation", "Angstrom", 350.0, [0, inf], "volume", 83 ["edge_separation", "Angstrom", 350.0, [0, inf], "volume", 84 84 "Mean separation of chained particles"], 85 ["string_thickness", "Angstrom", 2.5, [0, inf], "volume", 85 ["string_thickness", "Angstrom", 2.5, [0, inf], "volume", 86 86 "Thickness of the chain linkage"], 87 ["number_of_pearls", "none", 3, [0, inf], "volume", 87 ["number_of_pearls", "none", 3, [0, inf], "volume", 88 88 "Mean number of pearls in each necklace"], 89 ["sld", "Angstrom^2", 1.0, [-inf, inf], "", 89 ["sld", "Angstrom^2", 1.0, [-inf, inf], "", 90 90 "Scattering length density of the chained spheres"], 91 ["string_sld", "Angstrom^2", 1.0, [-inf, inf], "", 91 ["string_sld", "Angstrom^2", 1.0, [-inf, inf], "", 92 92 "Scattering length density of the chain linkage"], 93 ["solvent_sld", "Angstrom^2", 6.3, [-inf, inf], "", 93 ["solvent_sld", "Angstrom^2", 6.3, [-inf, inf], "", 94 94 "Scattering length density of the solvent"], 95 95 ] 96 96 97 97 source = ["lib/Si.c", "pearl_necklace.c"] 98 single = False # use double precision unless told otherwise 99 100 def volume(radius, edge_separation, string_thickness, number_of_pearls): 101 """ 102 Calculates the total particle volume of the necklace. 103 Redundant with form_volume. 104 """ 105 number_of_strings = number_of_pearls - 1.0 106 string_vol = edge_separation * pi * pow((string_thickness / 2.0), 2.0) 107 pearl_vol = 4.0 /3.0 * pi * pow(radius, 3.0) 108 total_vol = number_of_strings * string_vol 109 total_vol += number_of_pearls * pearl_vol 110 return total_vol 111 112 def ER(radius, edge_separation, string_thickness, number_of_pearls): 113 """ 114 Calculation for effective radius. 115 """ 116 tot_vol = volume(radius, edge_separation, string_thickness, number_of_pearls) 117 rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333) 118 return rad_out 98 119 99 120 # parameters for demo … … 110 131 # names and the target sasview model name. 111 132 oldname = 'PearlNecklaceModel' 112 oldpars = dict(scale='scale', background='background',radius='radius',133 oldpars = dict(scale='scale', background='background', radius='radius', 113 134 number_of_pearls='num_pearls', solvent_sld='sld_solv', 114 135 string_thickness='thick_string', sld='sld_pearl', 115 136 string_sld='sld_string', edge_separation='edge_separation') 137 138 tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]] -
sasmodels/models/polymer_excl_volume.py
r513efc5 r168052c 113 113 category = "shape-independent" 114 114 115 # pylint: disable=bad-whitespace, line-too-long 115 116 # ["name", "units", default, [lower, upper], "type", "description"], 116 117 parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"], 117 118 ["porod_exp", "", 3.0, [-inf, inf], "", "Porod exponent"], 118 ] 119 ] 120 # pylint: enable=bad-whitespace, line-too-long 119 121 120 122 121 def Iq(q, rg, porod_exp): 122 123 def Iq(q, 124 rg=60.0, 125 porod_exp=3.0): 123 126 """ 124 127 :param q: Input q-value (float or [float, float]) … … 131 134 o2nu = 1.0/(2.0*nu) 132 135 133 intensity = ((1.0/(nu*power(u, o2nu))) * (gamma(o2nu)*gammainc(o2nu, u) - 136 intensity = ((1.0/(nu*power(u, o2nu))) * 137 (gamma(o2nu)*gammainc(o2nu, u) - 134 138 1.0/power(u, o2nu) * gamma(porod_exp) * 135 139 gammainc(porod_exp, u))) * (q > 0) + 1.0*(q <= 0) … … 141 145 142 146 def Iqxy(qx, qy, *args): 143 iq = Iq(sqrt(qx**2 + qy**2), *args) 147 """ 148 :param qx: Input q_x-value 149 :param qy: Input q_y-value 150 :param args: Remaining arguments 151 :return: 2D-Intensity 152 """ 144 153 145 return iq154 return Iq(sqrt(qx**2 + qy**2), *args) 146 155 147 156 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values … … 157 166 porod_exp='m') 158 167 159 tests = [[{'rg': 10, 'porod_exp': 4.0}, 0.1, 0.723436675809], 168 tests = [ 169 # Accuracy tests based on content in test/polyexclvol_default_igor.txt 170 [{'rg': 60, 'porod_exp': 3.0}, 0.001, 0.998801], 171 [{'rg': 60, 'porod_exp': 3.0}, 0.105363, 0.0162751], 172 [{'rg': 60, 'porod_exp': 3.0}, 0.665075, 6.56261e-05], 160 173 161 [{'rg': 2.2, 'porod_exp': 22.0, 'background': 100.0}, 5.0, 100.0], 162 163 [{'rg': 1.1, 'porod_exp': 1, 'background': 10.0, 'scale': 1.25}, 164 20000., 10.0000712097] 165 ] 174 # Additional tests with larger range of parameters 175 [{'rg': 10, 'porod_exp': 4.0}, 0.1, 0.723436675809], 176 [{'rg': 2.2, 'porod_exp': 22.0, 'background': 100.0}, 5.0, 100.0], 177 [{'rg': 1.1, 'porod_exp': 1, 'background': 10.0, 'scale': 1.25}, 178 20000., 10.0000712097] 179 ] -
sasmodels/models/power_law.py
reb69cce r841753c 12 12 I(q) = \text{scale} \cdot q^{-\text{power}} + \text{background} 13 13 14 Note the minus sign in front of the exponent. The exponent *power* 14 Note the minus sign in front of the exponent. The exponent *power* 15 15 should therefore be entered as a **positive** number for fitting. 16 16 17 Also note that unlike many other models, *scale* in this model 18 is NOT explicitly related to a volume fraction. Be careful if 17 Also note that unlike many other models, *scale* in this model 18 is NOT explicitly related to a volume fraction. Be careful if 19 19 combining this model with other models. 20 20 … … 31 31 from numpy import inf, sqrt 32 32 33 name = 33 name = "power_law" 34 34 title = "Simple power law with a flat background" 35 35 36 description = """ \37 38 39 40 36 description = """ 37 Evaluates the function 38 I(q) = scale * q^(-power) + background 39 NB: enter power as a positive number! 40 """ 41 41 category = "shape-independent" 42 42 … … 45 45 46 46 # NB: Scale and Background are implicit parameters on every model 47 def Iq(q,power): 47 def Iq(q, power): 48 # pylint: disable=missing-docstring 48 49 inten = (q**-power) 49 50 return inten … … 51 52 52 53 def Iqxy(qx, qy, *args): 54 # pylint: disable=missing-docstring 53 55 return Iq(sqrt(qx ** 2 + qy ** 2), *args) 54 56 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values … … 64 66 65 67 tests = [ 66 [ {'scale': 1.0, 'power': 4.0, 'background' : 0.0}, [0.0106939, 0.469418], [7.64644e+07, 20.5949]] 67 ] 68 [{'scale': 1.0, 'power': 4.0, 'background' : 0.0}, 69 [0.0106939, 0.469418], [7.64644e+07, 20.5949]], 70 ] -
sasmodels/models/rpa.c
r82c299f r13ed84c 205 205 const double Kbb = 0.0; 206 206 const double Kcc = 0.0; 207 const double Kdd = 0.0;207 //const double Kdd = 0.0; 208 208 const double Zaa = Kaa - Kad - Kad; 209 209 const double Zab = Kab - Kad - Kbd; … … 278 278 const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ; 279 279 const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ; 280 const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;280 //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ; 281 281 const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ; 282 282 const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ; 283 const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;284 const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;283 //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ; 284 //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ; 285 285 const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ; 286 286 -
sasmodels/models/star_polymer.py
r55b283e8 r13ed84c 55 55 """ 56 56 category = "shape-independent" 57 57 single = False 58 # pylint: disable=bad-whitespace, line-too-long 58 59 # ["name", "units", default, [lower, upper], "type","description"], 59 60 parameters = [["radius2", "Ang", 100.0, [0.0, inf], "", "Ensemble radius of gyration squared of an arm"], 60 61 ["arms", "", 3, [1.0, 6.0], "", "Number of arms in the model"], 61 62 62 ] 63 # pylint: enable=bad-whitespace, line-too-long 63 64 64 65 source = ["star_polymer.c"] … … 75 76 tests = [[{'radius2': 2.0, 76 77 'arms': 3.3, 77 78 }, 0.5, 0.850646091108], 78 79 79 80 [{'radius2': 1.0, 80 81 'arms': 2.0, 81 82 'background': 1.8, 82 83 83 }, 1.0, 2.53575888234], 84 ] -
sasmodels/models/stickyhardsphere.py
rd18f8a8 r13ed84c 85 85 category = "structure-factor" 86 86 87 single = False 87 88 # ["name", "units", default, [lower, upper], "type","description"], 88 89 parameters = [ … … 182 183 demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05, 183 184 stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40) 185 # 186 tests = [ 187 [ {'scale': 1.0, 'background' : 0.0, 'effect_radius' : 50.0, 'perturb' : 0.05, 'stickiness' : 0.2, 'volfraction' : 0.1, 188 'effect_radius_pd' : 0}, [0.001], [1.09718]] 189 ] 184 190 191 -
sasmodels/models/surface_fractal.py
r9c461c7 r168052c 80 80 category = "shape-independent" 81 81 82 # pylint: disable=bad-whitespace, line-too-long 82 83 # ["name", "units", default, [lower, upper], "type","description"], 83 84 parameters = [["radius", "Ang", 10.0, [0, inf], "", … … 87 88 ["cutoff_length", "Ang", 500., [0.0, inf], "", 88 89 "Cut-off Length"], 89 90 90 ] 91 # pylint: enable=bad-whitespace, line-too-long 91 92 92 93 source = ["lib/sph_j1c.c", "lib/lanczos_gamma.c", "surface_fractal.c"] … … 100 101 cutoff_length='co_length') 101 102 102 tests = [[{'radius': 1.0, 'surface_dim': 1.0, 'cutoff_length': 10.0, 103 }, 0.332070182643, 1125.00321004], 103 tests = [ 104 # Accuracy tests based on content in test/utest_other_models.py 105 [{'radius': 10.0, 106 'surface_dim': 2.0, 107 'cutoff_length': 500.0, 108 }, 0.05, 301428.65916], 104 109 105 [{'radius': 3.5, 'surface_dim': 0.1, 'cutoff_length': 30.0, 106 'background': 0.01, 107 }, 5.0, 0.00999998891322], 110 # Additional tests with larger range of parameters 111 [{'radius': 1.0, 112 'surface_dim': 1.0, 113 'cutoff_length': 10.0, 114 }, 0.332070182643, 1125.00321004], 108 115 109 [{'radius': 3.0, 'surface_dim': 1.0, 'cutoff_length': 33.0, 110 'scale': 0.1, 111 }, 0.51, 2.50020147004], 112 ] 116 [{'radius': 3.5, 117 'surface_dim': 0.1, 118 'cutoff_length': 30.0, 119 'background': 0.01, 120 }, 5.0, 0.00999998891322], 121 122 [{'radius': 3.0, 123 'surface_dim': 1.0, 124 'cutoff_length': 33.0, 125 'scale': 0.1, 126 }, 0.51, 2.50020147004], 127 ] -
sasmodels/models/two_lorentzian.py
r513efc5 r168052c 54 54 category = "shape-independent" 55 55 56 # pylint: disable=bad-whitespace, line-too-long 56 57 # ["name", "units", default, [lower, upper], "type", "description"], 57 parameters = [["lorentz_scale_1", "", 10.0, [-inf, inf], "", 58 "First power law scale factor"], 59 ["lorentz_length_1", "Ang", 100.0, [-inf, inf], "", 60 "First Lorentzian screening length"], 61 ["lorentz_exp_1", "", 3.0, [-inf, inf], "", 62 "First exponent of power law"], 63 ["lorentz_scale_2", "", 1.0, [-inf, inf], "", 64 "Second scale factor for broad Lorentzian peak"], 65 ["lorentz_length_2", "Ang", 10.0, [-inf, inf], "", 66 "Second Lorentzian screening length"], 67 ["lorentz_exp_2", "", 2.0, [-inf, inf], "", 68 "Second exponent of power law"], 69 ] 58 parameters = [["lorentz_scale_1", "", 10.0, [-inf, inf], "", "First power law scale factor"], 59 ["lorentz_length_1", "Ang", 100.0, [-inf, inf], "", "First Lorentzian screening length"], 60 ["lorentz_exp_1", "", 3.0, [-inf, inf], "", "First exponent of power law"], 61 ["lorentz_scale_2", "", 1.0, [-inf, inf], "", "Second scale factor for broad Lorentzian peak"], 62 ["lorentz_length_2", "Ang", 10.0, [-inf, inf], "", "Second Lorentzian screening length"], 63 ["lorentz_exp_2", "", 2.0, [-inf, inf], "", "Second exponent of power law"], 64 ] 65 # pylint: enable=bad-whitespace, line-too-long 70 66 71 67 72 68 def Iq(q, 73 lorentz_scale_1 ,74 lorentz_length_1 ,75 lorentz_exp_1 ,76 lorentz_scale_2 ,77 lorentz_length_2 ,78 lorentz_exp_2 ):69 lorentz_scale_1=10.0, 70 lorentz_length_1=100.0, 71 lorentz_exp_1=3.0, 72 lorentz_scale_2=1.0, 73 lorentz_length_2=10.0, 74 lorentz_exp_2=2.0): 79 75 80 76 """ … … 88 84 :return: Calculated intensity 89 85 """ 90 86 # pylint: disable=bad-whitespace 91 87 intensity = lorentz_scale_1/(1.0 + 92 88 power(q*lorentz_length_1, lorentz_exp_1)) 93 89 intensity += lorentz_scale_2/(1.0 + 94 90 power(q*lorentz_length_2, lorentz_exp_2)) 91 # pylint: enable=bad-whitespace 95 92 96 93 return intensity … … 100 97 101 98 def Iqxy(qx, qy, *args): 102 iq = Iq(sqrt(qx**2 + qy**2), *args) 99 """ 100 :param qx: Input q_x-value 101 :param qy: Input q_y-value 102 :param args: Remaining arguments 103 :return: 2D-Intensity 104 """ 103 105 104 return iq106 return Iq(sqrt(qx**2 + qy**2), *args) 105 107 106 108 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values … … 108 110 109 111 demo = dict(scale=1, background=0.1, 110 lorentz_scale_1=10, lorentz_length_1=100.0, lorentz_exp_1=3.0, 111 lorentz_scale_2=1, lorentz_length_2=10, lorentz_exp_2=2.0) 112 lorentz_scale_1=10, 113 lorentz_length_1=100.0, 114 lorentz_exp_1=3.0, 115 lorentz_scale_2=1, 116 lorentz_length_2=10, 117 lorentz_exp_2=2.0) 112 118 113 119 oldname = "TwoLorentzianModel" 114 120 oldpars = dict(background='background', 115 lorentz_scale_1='scale_1', lorentz_scale_2='scale_2', 116 lorentz_length_1='length_1', lorentz_length_2='length_2', 117 lorentz_exp_1='exponent_1', lorentz_exp_2='exponent_2') 121 lorentz_scale_1='scale_1', 122 lorentz_scale_2='scale_2', 123 lorentz_length_1='length_1', 124 lorentz_length_2='length_2', 125 lorentz_exp_1='exponent_1', 126 lorentz_exp_2='exponent_2') 118 127 119 tests = [[{'lorentz_scale_1': 10.0, 120 'lorentz_length_1': 100.0, 121 'lorentz_exp_1': 3.0, 122 'lorentz_scale_2': 1.0, 123 'lorentz_length_2': 10.0, 124 'lorentz_exp_2': 2.0, 125 }, 0.000332070182643, 10.9996228107], 128 tests = [ 126 129 127 [{'lorentz_scale_1': 0.0, 128 'lorentz_length_1': 0.0, 129 'lorentz_exp_1': 0.0, 130 'lorentz_scale_2': 0.0, 131 'lorentz_length_2': 0.0, 132 'lorentz_exp_2': 0.0, 133 'background': 100.0 134 }, 5.0, 100.0], 130 # Accuracy tests based on content in test/utest_extra_models.py 131 [{'lorentz_scale_1': 10.0, 132 'lorentz_length_1': 100.0, 133 'lorentz_exp_1': 3.0, 134 'lorentz_scale_2': 1.0, 135 'lorentz_length_2': 10.0, 136 'lorentz_exp_2': 2.0, 137 'background': 0.1, 138 }, 0.001, 11.08991], 135 139 136 [{'lorentz_scale_1': 200.0, 137 'lorentz_length_1': 10.0, 138 'lorentz_exp_1': 0.1, 139 'lorentz_scale_2': 0.1, 140 'lorentz_length_2': 5.0, 141 'lorentz_exp_2': 2.0 142 }, 20000., 45.5659201896], 143 ] 140 [{'lorentz_scale_1': 10.0, 141 'lorentz_length_1': 100.0, 142 'lorentz_exp_1': 3.0, 143 'lorentz_scale_2': 1.0, 144 'lorentz_length_2': 10.0, 145 'lorentz_exp_2': 2.0, 146 'background': 0.1, 147 }, 0.150141, 0.410245], 148 149 [{'lorentz_scale_1': 10.0, 150 'lorentz_length_1': 100.0, 151 'lorentz_exp_1': 3.0, 152 'lorentz_scale_2': 1.0, 153 'lorentz_length_2': 10.0, 154 'lorentz_exp_2': 2.0, 155 'background': 0.1, 156 }, 0.442528, 0.148699], 157 158 # Additional tests with larger range of parameters 159 [{'lorentz_scale_1': 10.0, 160 'lorentz_length_1': 100.0, 161 'lorentz_exp_1': 3.0, 162 'lorentz_scale_2': 1.0, 163 'lorentz_length_2': 10.0, 164 'lorentz_exp_2': 2.0, 165 }, 0.000332070182643, 10.9996228107], 166 167 [{'lorentz_scale_1': 0.0, 168 'lorentz_length_1': 0.0, 169 'lorentz_exp_1': 0.0, 170 'lorentz_scale_2': 0.0, 171 'lorentz_length_2': 0.0, 172 'lorentz_exp_2': 0.0, 173 'background': 100.0 174 }, 5.0, 100.0], 175 176 [{'lorentz_scale_1': 200.0, 177 'lorentz_length_1': 10.0, 178 'lorentz_exp_1': 0.1, 179 'lorentz_scale_2': 0.1, 180 'lorentz_length_2': 5.0, 181 'lorentz_exp_2': 2.0 182 }, 20000., 45.5659201896], 183 ] -
sasmodels/resolution.py
r5258859 r5925e90 5 5 """ 6 6 from __future__ import division 7 8 from scipy.special import erf 9 from numpy import sqrt, log, log10 10 import numpy as np 7 11 8 12 __all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", … … 10 14 "pinhole_extend_q", "slit_extend_q", "bin_edges", 11 15 "interpolate", "linear_extrapolation", "geometric_extrapolation", 12 ] 13 14 from scipy.special import erf 15 from numpy import sqrt, log, log10 16 import numpy as np 16 ] 17 17 18 18 MINIMUM_RESOLUTION = 1e-8 … … 80 80 # default to Perfect1D if the pinhole geometry is not defined. 81 81 self.q, self.q_width = q, q_width 82 self.q_calc = pinhole_extend_q(q, q_width, nsigma=nsigma) \83 if q_calc is None else np.sort(q_calc)84 self.weight_matrix = pinhole_resolution( self.q_calc,85 82 self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma) 83 if q_calc is None else np.sort(q_calc)) 84 self.weight_matrix = pinhole_resolution( 85 self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 86 86 87 87 def apply(self, theory): … … 135 135 """ 136 136 #print("apply shapes", theory.shape, weight_matrix.shape) 137 Iq = np.dot(theory[None, :], weight_matrix)137 Iq = np.dot(theory[None, :], weight_matrix) 138 138 #print("result shape",Iq.shape) 139 139 return Iq.flatten() … … 153 153 # neither trapezoid nor Simpson's rule improved the accuracy. 154 154 edges = bin_edges(q_calc) 155 edges[edges <0.0] = 0.0 # clip edges below zero156 G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:])155 edges[edges < 0.0] = 0.0 # clip edges below zero 156 G = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 157 157 weights = G[1:] - G[:-1] 158 weights /= np.sum(weights, axis=0)[None, :]158 weights /= np.sum(weights, axis=0)[None, :] 159 159 return weights 160 160 … … 287 287 # The current algorithm is a midpoint rectangle rule. 288 288 q_edges = bin_edges(q_calc) # Note: requires q > 0 289 q_edges[q_edges <0.0] = 0.0 # clip edges below zero289 q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 290 290 weights = np.zeros((len(q), len(q_calc)), 'd') 291 291 … … 306 306 #print(qi - h, qi + h) 307 307 #print(in_x + abs_x) 308 weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)308 weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 309 309 else: 310 310 L = n_height 311 311 for k in range(-L, L+1): 312 weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w)313 weights[i, :] /= 2*L + 1312 weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w) 313 weights[i, :] /= 2*L + 1 314 314 315 315 return weights.T … … 358 358 log-scaled data. 359 359 """ 360 if len(x) < 2 or (np.diff(x) <0).any():360 if len(x) < 2 or (np.diff(x) < 0).any(): 361 361 raise ValueError("Expected bins to be an increasing set") 362 362 edges = np.hstack([ … … 373 373 """ 374 374 step = np.diff(q) 375 index = step >max_step375 index = step > max_step 376 376 if np.any(index): 377 377 inserts = [] 378 for q_i, step_i in zip(q[:-1][index],step[index]):378 for q_i, step_i in zip(q[:-1][index], step[index]): 379 379 n = np.ceil(step_i/max_step) 380 inserts.extend(q_i + np.arange(1, n)*(step_i/n))380 inserts.extend(q_i + np.arange(1, n)*(step_i/n)) 381 381 # Extend a couple of fringes beyond the end of the data 382 inserts.extend(q[-1] + np.arange(1, 8)*max_step)383 q_calc = np.sort(np.hstack((q, inserts)))382 inserts.extend(q[-1] + np.arange(1, 8)*max_step) 383 q_calc = np.sort(np.hstack((q, inserts))) 384 384 else: 385 385 q_calc = q … … 399 399 if q_min < q[0]: 400 400 if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 401 n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] >q[0] else 15401 n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 402 402 q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 403 403 else: 404 404 q_low = [] 405 405 if q_max > q[-1]: 406 n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] >q[-2] else 15406 n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15 407 407 q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 408 408 else: … … 452 452 if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 453 453 n_low = log_delta_q * (log(q[0])-log(q_min)) 454 q_low 454 q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 455 455 else: 456 456 q_low = [] … … 470 470 471 471 def eval_form(q, form, pars): 472 """ 473 Return the SAS model evaluated at *q*. 474 475 *form* is the SAS model returned from :fun:`core.load_model`. 476 477 *pars* are the parameter values to use when evaluating. 478 """ 472 479 from sasmodels import core 473 480 kernel = core.make_kernel(form, [q]) … … 478 485 479 486 def gaussian(q, q0, dq): 487 """ 488 Return the Gaussian resolution function. 489 490 *q0* is the center, *dq* is the width and *q* are the points to evaluate. 491 """ 480 492 from numpy import exp, pi 481 493 return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) … … 489 501 that make it slow to evaluate but give it good accuracy. 490 502 """ 491 from scipy.integrate import romberg , dblquad503 from scipy.integrate import romberg 492 504 493 505 if any(k not in form.info['defaults'] for k in pars.keys()): … … 520 532 result[i] = r/(2*h) 521 533 else: 522 w_grid = np.linspace(0, w, 21)[None, :]523 h_grid = np.linspace(-h, h, 23)[:, None]534 w_grid = np.linspace(0, w, 21)[None, :] 535 h_grid = np.linspace(-h, h, 23)[:, None] 524 536 u = sqrt((qi+h_grid)**2 + w_grid**2) 525 537 Iu = np.interp(u, q_calc, Iq) 526 538 #print(np.trapz(Iu, w_grid, axis=1)) 527 Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0])539 Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0]) 528 540 result[i] = Is / (2*h*w) 529 """ 530 r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, 531 args=(qi,)) 532 result[i] = r/(w*2*h) 533 """ 541 # from scipy.integrate import dblquad 542 # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, 543 # args=(qi,)) 544 # result[i] = r/(w*2*h) 534 545 535 546 # r should be [float, ...], but it is [array([float]), array([float]),...] … … 553 564 554 565 _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 555 r = [romberg(_fn, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi, args=(qi, dqi),556 divmax=100, vec_func=True, tol=0, rtol=1e-8)557 for qi, dqi in zip(q,q_width)]566 r = [romberg(_fn, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi, 567 args=(qi, dqi), divmax=100, vec_func=True, tol=0, rtol=1e-8) 568 for qi, dqi in zip(q, q_width)] 558 569 return np.asarray(r).flatten() 559 570 560 571 561 572 class ResolutionTest(unittest.TestCase): 573 """ 574 Test the resolution calculations. 575 """ 562 576 563 577 def setUp(self): … … 595 609 theory = self.Iq(resolution.q_calc) 596 610 output = resolution.apply(theory) 597 answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 598 5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ] 611 answer = [ 612 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 613 5.5555, 4.5584, 3.5606, 2.5623, 2.0000, 614 ] 599 615 np.testing.assert_allclose(output, answer, atol=1e-4) 600 616 … … 608 624 theory = 1000*self.Iq(resolution.q_calc**4) 609 625 output = resolution.apply(theory) 610 answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 611 5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ] 626 answer = [ 627 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 628 5.40363, 4.40655, 3.40880, 2.41058, 2.00000, 629 ] 612 630 np.testing.assert_allclose(output, answer, atol=1e-4) 613 631 … … 620 638 theory = self.Iq(resolution.q_calc) 621 639 output = resolution.apply(theory) 622 answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 640 answer = [ 641 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 642 ] 623 643 np.testing.assert_allclose(output, answer, atol=1e-4) 624 644 … … 632 652 theory = self.Iq(resolution.q_calc) 633 653 output = resolution.apply(theory) 634 answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 654 answer = [ 655 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 656 ] 635 657 np.testing.assert_allclose(output, answer, atol=1e-4) 636 658 … … 652 674 theory = 12.0-1000.0*resolution.q_calc 653 675 output = resolution.apply(theory) 654 answer = [ 10.44785079, 9.84991299, 8.98101708, 655 7.99906585, 6.99998311, 6.00001689, 656 5.00093415, 4.01898292, 3.15008701, 2.55214921] 676 answer = [ 677 10.44785079, 9.84991299, 8.98101708, 678 7.99906585, 6.99998311, 6.00001689, 679 5.00093415, 4.01898292, 3.15008701, 2.55214921, 680 ] 657 681 np.testing.assert_allclose(output, answer, atol=1e-8) 658 682 659 683 660 684 class IgorComparisonTest(unittest.TestCase): 685 """ 686 Test resolution calculations against those returned by Igor. 687 """ 661 688 662 689 def setUp(self): … … 666 693 self.model = core.load_model(sphere, dtype='double') 667 694 668 def Iq_sphere(self, pars, resolution):695 def _eval_sphere(self, pars, resolution): 669 696 from sasmodels import core 670 697 kernel = core.make_kernel(self.model, [resolution.q_calc]) … … 674 701 return result 675 702 676 def compare(self, q, output, answer, tolerance):703 def _compare(self, q, output, answer, tolerance): 677 704 #err = (output - answer)/answer 678 705 #idx = abs(err) >= tolerance … … 691 718 q, width, answer, _ = data 692 719 resolution = Perfect1D(q) 693 output = self. Iq_sphere(pars, resolution)694 self. compare(q, output, answer, 1e-6)720 output = self._eval_sphere(pars, resolution) 721 self._compare(q, output, answer, 1e-6) 695 722 696 723 def test_pinhole(self): … … 704 731 q, q_width, answer = data 705 732 resolution = Pinhole1D(q, q_width) 706 output = self. Iq_sphere(pars, resolution)733 output = self._eval_sphere(pars, resolution) 707 734 # TODO: relative error should be lower 708 self. compare(q, output, answer, 3e-4)735 self._compare(q, output, answer, 3e-4) 709 736 710 737 def test_pinhole_romberg(self): … … 727 754 q_calc, tol = None, 0.01 728 755 resolution = Pinhole1D(q, q_width, q_calc=q_calc) 729 output = self. Iq_sphere(pars, resolution)756 output = self._eval_sphere(pars, resolution) 730 757 if 0: # debug plot 731 758 import matplotlib.pyplot as plt 732 759 resolution = Perfect1D(q) 733 source = self. Iq_sphere(pars, resolution)760 source = self._eval_sphere(pars, resolution) 734 761 plt.loglog(q, source, '.') 735 762 plt.loglog(q, answer, '-', hold=True) 736 763 plt.loglog(q, output, '-', hold=True) 737 764 plt.show() 738 self. compare(q, output, answer, tol)765 self._compare(q, output, answer, tol) 739 766 740 767 def test_slit(self): … … 748 775 q, delta_qv, _, answer = data 749 776 resolution = Slit1D(q, width=delta_qv, height=0) 750 output = self. Iq_sphere(pars, resolution)777 output = self._eval_sphere(pars, resolution) 751 778 # TODO: eliminate Igor test since it is too inaccurate to be useful. 752 779 # This means we can eliminate the test data as well, and instead 753 780 # use a generated q vector. 754 self. compare(q, output, answer, 0.5)781 self._compare(q, output, answer, 0.5) 755 782 756 783 def test_slit_romberg(self): … … 768 795 delta_qv[0], 0.) 769 796 resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc) 770 output = self. Iq_sphere(pars, resolution)797 output = self._eval_sphere(pars, resolution) 771 798 # TODO: relative error should be lower 772 self. compare(q, output, answer, 0.025)799 self._compare(q, output, answer, 0.025) 773 800 774 801 def test_ellipsoid(self): … … 783 810 } 784 811 form = load_model('ellipsoid', dtype='double') 785 q = np.logspace(log10(4e-5), log10(2.5e-2), 68)812 q = np.logspace(log10(4e-5), log10(2.5e-2), 68) 786 813 width, height = 0.117, 0. 787 814 resolution = Slit1D(q, width=width, height=height) … … 790 817 # TODO: 10% is too much error; use better algorithm 791 818 #print(np.max(abs(answer-output)/answer)) 792 self. compare(q, output, answer, 0.1)819 self._compare(q, output, answer, 0.1) 793 820 794 821 #TODO: can sas q spacing be too sparse for the resolution calculation? … … 804 831 q, q_width, answer = data[:, ::20] # Take every nth point 805 832 resolution = Pinhole1D(q, q_width) 806 output = self. Iq_sphere(pars, resolution)807 self. compare(q, output, answer, 1e-6)833 output = self._eval_sphere(pars, resolution) 834 self._compare(q, output, answer, 1e-6) 808 835 809 836 … … 1058 1085 1059 1086 def demo_pinhole_1d(): 1087 """ 1088 Show example of pinhole smearing. 1089 """ 1060 1090 q = np.logspace(-4, np.log10(0.2), 400) 1061 1091 q_width = 0.1*q … … 1064 1094 1065 1095 def demo_slit_1d(): 1096 """ 1097 Show example of slit smearing. 1098 """ 1066 1099 q = np.logspace(-4, np.log10(0.2), 100) 1067 1100 w = h = 0. … … 1071 1104 #h = 0.0277790 1072 1105 resolution = Slit1D(q, w, h) 1073 _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))1106 _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h)) 1074 1107 1075 1108 def demo(): 1109 """ 1110 Run the resolution demos. 1111 """ 1076 1112 import matplotlib.pyplot as plt 1077 1113 plt.subplot(121) -
sasmodels/resolution2d.py
r9404dd3 r823e620 2 2 #This software was developed by the University of Tennessee as part of the 3 3 #Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 4 #project funded by the US National Science Foundation. 4 #project funded by the US National Science Foundation. 5 5 #See the license text in license.txt 6 6 """ … … 19 19 ## Defaults 20 20 NR = {'xhigh':10, 'high':5, 'med':5, 'low':3} 21 NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4}21 NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4} 22 22 23 23 class Pinhole2D(Resolution): … … 25 25 Gaussian Q smearing class for SAS 2d data 26 26 """ 27 27 28 28 def __init__(self, data=None, index=None, 29 29 nsigma=NSIGMA, accuracy='Low', coords='polar'): 30 30 """ 31 31 Assumption: equally spaced bins in dq_r, dq_phi space. 32 32 33 33 :param data: 2d data used to set the smearing parameters 34 34 :param index: 1d array with len(data) to define the range … … 42 42 ## number of bins in r axis for over-sampling 43 43 self.nr = NR[accuracy.lower()] 44 ## number of bins in phi axis for over-sampling 44 ## number of bins in phi axis for over-sampling 45 45 self.nphi = NPHI[accuracy.lower()] 46 46 ## maximum nsigmas 47 self.nsigma = nsigma47 self.nsigma = nsigma 48 48 self.coords = coords 49 49 self._init_data(data, index) … … 85 85 def _calc_res(self): 86 86 """ 87 Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, 87 Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, 88 88 then find smeared intensity 89 """ 89 """ 90 90 nr, nphi = self.nr, self.nphi 91 91 # Total number of bins = # of bins … … 143 143 if self.coords == 'polar': 144 144 q_r = sqrt(qx**2 + qy**2) 145 qx_res = ( (dqx*cos(dphi) + q_r) * cos(-q_phi) +146 147 qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) +148 145 qx_res = ((dqx*cos(dphi) + q_r) * cos(-q_phi) 146 + dqy*sin(dphi) * sin(-q_phi)) 147 qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) 148 + dqy*sin(dphi) * cos(-q_phi)) 149 149 else: 150 qx_res = qx + 151 qy_res = qy + 150 qx_res = qx + dqx*cos(dphi) 151 qy_res = qy + dqy*sin(dphi) 152 152 153 153 … … 162 162 theory = np.reshape(theory, (nbins, nq)) 163 163 ## Averaging with Gaussian weighting: normalization included. 164 value = np.average(theory, axis=0, weights=self.q_calc_weights)164 value = np.average(theory, axis=0, weights=self.q_calc_weights) 165 165 ## Return the smeared values in the range of self.index 166 166 return value 167 167 else: 168 168 return theory 169 170 """171 if __name__ == '__main__':172 ## Test w/ 2D linear function173 x = 0.001*np.arange(1, 11)174 dx = np.ones(len(x))*0.0003175 y = 0.001*np.arange(1, 11)176 dy = np.ones(len(x))*0.001177 z = np.ones(10)178 dz = sqrt(z)179 180 from sas.dataloader import Data2D181 #for i in range(10): print(i, 0.001 + i*0.008/9.0)182 #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) )))183 out = Data2D()184 out.data = z185 out.qx_data = x186 out.qy_data = y187 out.dqx_data = dx188 out.dqy_data = dy189 out.q_data = sqrt(dx * dx + dy * dy)190 index = np.ones(len(x), dtype = bool)191 out.mask = index192 from sas.models.LineModel import LineModel193 model = LineModel()194 model.setParam("A", 0)195 196 smear = Smearer2D(out, model, index)197 #smear.set_accuracy('Xhigh')198 value = smear.get_value()199 ## All data are ones, so the smeared should also be ones.200 print("Data length =", len(value))201 print(" 2D linear function, I = 0 + 1*qy")202 text = " Gaussian weighted averaging on a 2D linear function will "203 text += "provides the results same as without the averaging."204 print(text)205 print("qx_data", "qy_data", "I_nonsmear", "I_smeared")206 for ind in range(len(value)):207 print(x[ind], y[ind], model.evalDistribution([x, y])[ind], value[ind])208 209 210 if __name__ == '__main__':211 ## Another Test w/ constant function212 x = 0.001*np.arange(1,11)213 dx = np.ones(len(x))*0.001214 y = 0.001*np.arange(1,11)215 dy = np.ones(len(x))*0.001216 z = np.ones(10)217 dz = sqrt(z)218 219 from DataLoader import Data2D220 #for i in range(10): print(i, 0.001 + i*0.008/9.0)221 #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) )))222 out = Data2D()223 out.data = z224 out.qx_data = x225 out.qy_data = y226 out.dqx_data = dx227 out.dqy_data = dy228 index = np.ones(len(x), dtype = bool)229 out.mask = index230 from sas.models.Constant import Constant231 model = Constant()232 233 value = Smearer2D(out,model,index).get_value()234 ## All data are ones, so the smeared values should also be ones.235 print("Data length =",len(value), ", Data=",value)236 """ -
sasmodels/sasview_model.py
r9404dd3 reafc9fa 283 283 """ 284 284 q_vectors = [np.asarray(q) for q in args] 285 fn = self._model( self._model.make_input(q_vectors))285 fn = self._model(q_vectors) 286 286 pars = [self.params[v] for v in fn.fixed_pars] 287 287 pd_pars = [self._get_weights(p) for p in fn.pd_pars] -
sasmodels/sesans.py
r3c56da87 r190fc2b 12 12 import numpy as np 13 13 from numpy import pi, exp 14 15 14 from scipy.special import jv as besselj 16 15 17 def make_q(q_zmax, Rmax): 16 def make_q(q_max, Rmax): 17 r""" 18 Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 19 to $q_max$. 20 """ 18 21 q_min = dq = 0.1 * 2*pi / Rmax 19 #q_min = 0.00003 20 return np.arange(q_min, q_zmax, dq) 21 22 # TODO: dead code; for now the call to the hankel transform happens in BumpsModel 23 class SesansCalculator: 24 def __init__(self, kernel, q_zmax, Rmax, SElength, wavelength, thickness): 25 self._set_kernel(kernel, q_zmax, Rmax) 26 self.SElength = SElength 27 self.wavelength = wavelength 28 self.thickness = thickness 29 30 def _set_kernel(self, kernel, q_zmax, Rmax): 31 kernel_input = kernel.make_input([make_q(q_zmax, Rmax)]) 32 self.sans_calculator = kernel(kernel_input) 33 34 def __call__(self, pars, pd_pars, cutoff=1e-5): 35 Iq = self.sans_calculator(pars, pd_pars, cutoff) 36 P = hankel(self.SElength, self.wavelength, self.thickness, self.q, Iq) 37 self.Iq = Iq 38 return P 22 return np.arange(q_min, q_max, dq) 39 23 40 24 def hankel(SElength, wavelength, thickness, q, Iq): 41 """25 r""" 42 26 Compute the expected SESANS polarization for a given SANS pattern. 43 27 44 Uses the hankel transform followed by the exponential. The values 45 for zz (or spin echo length, or delta), wavelength and sample thickness46 information should come from the dataset. *q* should be chosen such47 that the oscillations in *I(q)* are well sampled (e.g., 5*2*pi/d_max).28 Uses the hankel transform followed by the exponential. The values for *zz* 29 (or spin echo length, or delta), wavelength and sample thickness should 30 come from the dataset. $q$ should be chosen such that the oscillations 31 in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 48 32 49 *SElength* [A] is the set of z points at which to compute the hankel transform 33 *SElength* [A] is the set of $z$ points at which to compute the 34 Hankel transform 50 35 51 36 *wavelength* [m] is the wavelength of each individual point *zz* … … 53 38 *thickness* [cm] is the sample thickness. 54 39 55 *q* [A ^{-1}] is the set of q points at which the model has been computed.56 These should be equally spaced.40 *q* [A$^{-1}$] is the set of $q$ points at which the model has been 41 computed. These should be equally spaced. 57 42 58 *I* [cm ^{-1}] is the value of the SANS model at *q*43 *I* [cm$^{-1}$] is the value of the SANS model at *q* 59 44 """ 60 G = np.zeros(len(SElength), 'd') 61 for i in range(len(SElength)): 62 integr = besselj(0,q*SElength[i])*Iq*q 63 G[i] = np.sum(integr) 64 dq=(q[1]-q[0])*1e10 # [m^-1] step size in q, needed for integration 65 G *= dq*1e10*2*pi # integr step, conver q into [m**-1] and 2 pi circle integr 45 G = np.zeros_like(SElength, 'd') 46 for i, SElength_i in enumerate(SElength): 47 integral = besselj(0, q*SElength_i)*Iq*q 48 G[i] = np.sum(integral) 49 50 # [m^-1] step size in q, needed for integration 51 dq = (q[1]-q[0])*1e10 52 53 # integration step, convert q into [m**-1] and 2 pi circle integration 54 G *= dq*1e10*2*pi 55 66 56 P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 67 57 -
sasmodels/weights.py
r3c56da87 r5c962df 22 22 23 23 def get_pars(self): 24 """ 25 Return the parameters to the disperser as a dictionary. 26 """ 24 27 pars = {'type': self.type} 25 28 pars.update(self.__dict__) … … 28 31 # pylint: disable=no-self-use 29 32 def set_weights(self, values, weights): 33 """ 34 Set the weights on the disperser if it is :class:`ArrayDispersion`. 35 """ 30 36 raise RuntimeError("set_weights is only available for ArrayDispersion") 31 37 … … 36 42 *center* is the center of the distribution 37 43 38 *lb*, *ub* are the min and max allowed values44 *lb*, *ub* are the min and max allowed values 39 45 40 46 *relative* is True if the distribution width is proportional to the … … 63 69 64 70 class GaussianDispersion(Dispersion): 71 r""" 72 Gaussian dispersion, with 1-\ $\sigma$ width. 73 74 .. math:: 75 76 w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right) 77 """ 65 78 type = "gaussian" 66 79 default = dict(npts=35, width=0, nsigmas=3) … … 72 85 73 86 class RectangleDispersion(Dispersion): 87 r""" 88 Uniform dispersion, with width $\sqrt{3}\sigma$. 89 90 .. math:: 91 92 w = 1 93 """ 74 94 type = "rectangle" 75 95 default = dict(npts=35, width=0, nsigmas=1.70325) … … 81 101 82 102 class LogNormalDispersion(Dispersion): 103 r""" 104 log Gaussian dispersion, with 1-\ $\sigma$ width. 105 106 .. math:: 107 108 w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma} 109 """ 83 110 type = "lognormal" 84 111 default = dict(npts=80, width=0, nsigmas=8) 85 112 def _weights(self, center, sigma, lb, ub): 86 x = self._linspace(center, sigma, max(lb, 1e-8), max(ub,1e-8))87 px = np.exp(-0.5*(np.log(x)-center)**2 )/sigma**2/(x*sigma)113 x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) 114 px = np.exp(-0.5*(np.log(x)-center)**2/sigma**2)/(x*sigma) 88 115 return x, px 89 116 90 117 91 118 class SchulzDispersion(Dispersion): 119 r""" 120 Schultz dispersion, with 1-\ $\sigma$ width. 121 122 .. math:: 123 124 w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)} 125 126 where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$. 127 128 This is evaluated using logarithms as 129 130 .. math:: 131 132 w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right) 133 """ 92 134 type = "schulz" 93 135 default = dict(npts=80, width=0, nsigmas=8) 94 136 def _weights(self, center, sigma, lb, ub): 95 x = self._linspace(center, sigma, max(lb, 1e-8), max(ub,1e-8))96 R = x/center137 x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) 138 R = x/center 97 139 z = (center/sigma)**2 98 140 arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z) … … 102 144 103 145 class ArrayDispersion(Dispersion): 146 r""" 147 Empirical dispersion curve. 148 149 Use :meth:`set_weights` to set $w = f(x)$. 150 """ 104 151 type = "array" 105 152 default = dict(npts=35, width=0, nsigmas=1) … … 110 157 111 158 def set_weights(self, values, weights): 159 """ 160 Set the weights for the given x values. 161 """ 112 162 self.values = np.ascontiguousarray(values, 'd') 113 163 self.weights = np.ascontiguousarray(weights, 'd') … … 117 167 # TODO: interpolate into the array dispersion using npts 118 168 x = center + self.values*sigma 119 idx = (x >=lb)&(x<=ub)169 idx = (x >= lb) & (x <= ub) 120 170 x = x[idx] 121 171 px = self.weights[idx] … … 124 174 125 175 # dispersion name -> disperser lookup table. 126 models = dict((d.type, d) for d in (176 models = dict((d.type, d) for d in ( 127 177 GaussianDispersion, RectangleDispersion, 128 178 ArrayDispersion, SchulzDispersion, LogNormalDispersion … … 149 199 of the parameter, and false if it is an absolute width. 150 200 151 Returns *(value, weight)*, where *value* and *weight* are vectors.201 Returns *(value, weight)*, where *value* and *weight* are vectors. 152 202 """ 153 203 cls = models[disperser] 154 204 obj = cls(n, width, nsigmas) 155 v, w = obj.get_weights(value, limits[0], limits[1], relative)156 return v, w205 v, w = obj.get_weights(value, limits[0], limits[1], relative) 206 return v, w 157 207 158 208 # Hack to allow sasview dispersion objects to interoperate with sasmodels 159 dispersers = dict((v.__name__, k) for k,v in models.items())209 dispersers = dict((v.__name__, k) for k, v in models.items()) 160 210 dispersers['DispersionModel'] = RectangleDispersion.type 161 211
Note: See TracChangeset
for help on using the changeset viewer.