Changes in / [204fd9b:5b0335b] in sasmodels
- Files:
-
- 6 added
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/developer/index.rst
r56fc97a rb85be2d 3 3 *************************** 4 4 5 .. toctree:: 6 :numbered: 4 7 :maxdepth: 4 5 8 9 calculator.rst -
sasmodels/bumps_model.py
rea75043 rea75043 81 81 from bumps.names import Parameter 82 82 83 pars = {} 83 pars = {} # floating point parameters 84 pd_types = {} # distribution names 84 85 for p in model_info['parameters']: 85 86 value = kwargs.pop(p.name, p.default) 86 87 pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 87 for name in model_info['partype']['pd-2d']: 88 for xpart, xdefault, xlimits in [ 89 ('_pd', 0., pars[name].limits), 90 ('_pd_n', 35., (0, 1000)), 91 ('_pd_nsigma', 3., (0, 10)), 92 ]: 93 xname = name + xpart 94 xvalue = kwargs.pop(xname, xdefault) 95 pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 96 97 pd_types = {} 98 for name in model_info['partype']['pd-2d']: 99 xname = name + '_pd_type' 100 xvalue = kwargs.pop(xname, 'gaussian') 101 pd_types[xname] = xvalue 88 if p.polydisperse: 89 for part, default, limits in [ 90 ('_pd', 0., pars[p.name].limits), 91 ('_pd_n', 35., (0, 1000)), 92 ('_pd_nsigma', 3., (0, 10)), 93 ]: 94 name = p.name + part 95 value = kwargs.pop(name, default) 96 pars[name] = Parameter.default(value, name=name, limits=limits) 97 pd_types[p.name+'_pd_type'] = kwargs.pop(name, 'gaussian') 102 98 103 99 if kwargs: # args not corresponding to parameters -
sasmodels/compare.py
rea75043 rea75043 306 306 Format the parameter list for printing. 307 307 """ 308 if is2d:309 exclude = lambda n: False310 else:311 partype = model_info['partype']312 par1d = set(partype['fixed-1d']+partype['pd-1d'])313 exclude = lambda n: n not in par1d314 308 lines = [] 315 for p in model_info['parameters']:316 if exclude(p.name): continue309 parameters = model_info['parameters'] 310 for p in parameters.user_parameters(pars, is2d): 317 311 fields = dict( 318 value=pars.get(p. name, p.default),319 pd=pars.get(p. name+"_pd", 0.),320 n=int(pars.get(p. name+"_pd_n", 0)),321 nsigma=pars.get(p. name+"_pd_nsgima", 3.),322 type=pars.get(p. name+"_pd_type", 'gaussian'))312 value=pars.get(p.id, p.default), 313 pd=pars.get(p.id+"_pd", 0.), 314 n=int(pars.get(p.id+"_pd_n", 0)), 315 nsigma=pars.get(p.id+"_pd_nsgima", 3.), 316 type=pars.get(p.id+"_pd_type", 'gaussian')) 323 317 lines.append(_format_par(p.name, **fields)) 324 318 return "\n".join(lines) … … 454 448 """ 455 449 # initialize the code so time is more accurate 456 value = calculator(**suppress_pd(pars)) 450 if Nevals > 1: 451 value = calculator(**suppress_pd(pars)) 457 452 toc = tic() 458 453 for _ in range(max(Nevals, 1)): # make sure there is at least one eval … … 661 656 """ 662 657 # Get the default values for the parameters 663 pars = dict((p.name, p.default) for p in model_info['parameters']) 664 665 # Fill in default values for the polydispersity parameters 666 for p in model_info['parameters']: 667 if p.type in ('volume', 'orientation'): 668 pars[p.name+'_pd'] = 0.0 669 pars[p.name+'_pd_n'] = 0 670 pars[p.name+'_pd_nsigma'] = 3.0 671 pars[p.name+'_pd_type'] = "gaussian" 658 pars = {} 659 for p in model_info['parameters'].call_parameters: 660 parts = [('', p.default)] 661 if p.polydisperse: 662 parts.append(('_pd', 0.0)) 663 parts.append(('_pd_n', 0)) 664 parts.append(('_pd_nsigma', 3.0)) 665 parts.append(('_pd_type', "gaussian")) 666 for ext, val in parts: 667 if p.length > 1: 668 dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length)) 669 else: 670 pars[p.id+ext] = val 672 671 673 672 # Plug in values given in demo … … 774 773 775 774 if len(engines) == 0: 776 engines.extend(['single', ' sasview'])775 engines.extend(['single', 'double']) 777 776 elif len(engines) == 1: 778 if engines[0][0] != ' sasview':779 engines.append(' sasview')777 if engines[0][0] != 'double': 778 engines.append('double') 780 779 else: 781 780 engines.append('single') … … 880 879 pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 881 880 if not opts['is2d']: 882 active = [base + ext 883 for base in model_info['partype']['pd-1d'] 884 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 885 active.extend(model_info['partype']['fixed-1d']) 886 for k in active: 887 v = pars[k] 888 v.range(*parameter_range(k, v.value)) 881 for p in model_info['parameters'].type['1d']: 882 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 883 k = p.name+ext 884 v = pars.get(k, None) 885 if v is not None: 886 v.range(*parameter_range(k, v.value)) 889 887 else: 890 888 for k, v in pars.items(): -
sasmodels/convert.py
r2a3e1f5 r2a3e1f5 73 73 new model definition end with sld. 74 74 """ 75 print "pars",pars 75 76 return dict((p, (v*1e-6 if p.startswith('sld') or p.endswith('sld') 76 77 else v*1e15 if 'ndensity' in p -
sasmodels/core.py
rb7172bb rb7172bb 6 6 "build_model", "make_kernel", "call_kernel", "call_ER_VR", 7 7 ] 8 9 8 10 9 from os.path import basename, dirname, join as joinpath, splitext … … 30 29 31 30 try: 31 # Python 3.5 and up 32 from importlib.util import spec_from_file_location, module_from_spec 33 def load_module(fullname, path): 34 spec = spec_from_file_location(fullname, path) 35 module = module_from_spec(spec) 36 spec.loader.exec_module(module) 37 return module 38 except ImportError: 39 # CRUFT: python 2 40 import imp 41 def load_module(fullname, path): 42 module = imp.load_source(fullname, path) 43 return module 44 45 try: 32 46 np.meshgrid([]) 33 47 meshgrid = np.meshgrid … … 40 54 return [np.asarray(v) for v in args] 41 55 42 43 56 def list_models(): 44 57 """ … … 63 76 """ 64 77 return build_model(load_model_info(model_name), **kw) 65 66 def load_model_info_from_path(path):67 # Pull off the last .ext if it exists; there may be others68 name = basename(splitext(path)[0])69 70 # Not cleaning name since don't need to be able to reload this71 # model later72 # Should probably turf the model from sys.modules after we are done...73 74 # Placing the model in the 'sasmodels.custom' name space, even75 # though it doesn't actually exist. imp.load_source doesn't seem76 # to care.77 kernel_module = imp.load_source('sasmodels.custom.'+name, path)78 79 # Now that we have the module, we can load it as usual80 model_info = generate.make_model_info(kernel_module)81 return model_info82 78 83 79 def load_model_info(model_name): … … 102 98 return product.make_product_info(P_info, Q_info) 103 99 100 kernel_module = load_kernel_module(model_name) 101 return generate.make_model_info(kernel_module) 102 103 104 def load_kernel_module(model_name): 105 if model_name.endswith('.py'): 106 path = model_name 107 # Pull off the last .ext if it exists; there may be others 108 name = basename(splitext(path)[0]) 109 # Placing the model in the 'sasmodels.custom' name space. 110 from sasmodels import custom 111 kernel_module = load_module('sasmodels.custom.'+name, path) 112 else: 113 from sasmodels import models 114 __import__('sasmodels.models.'+model_name) 115 kernel_module = getattr(models, model_name, None) 104 116 #import sys; print "\n".join(sys.path) 105 117 __import__('sasmodels.models.'+model_name) 106 118 kernel_module = getattr(models, model_name, None) 107 return generate.make_model_info(kernel_module)119 return kernel_module 108 120 109 121 … … 179 191 180 192 181 def make_kernel(model, q_vectors): 182 """ 183 Return a computation kernel from the model definition and the q input. 184 """ 185 return model(q_vectors) 186 187 def get_weights(model_info, pars, name): 193 def get_weights(parameter, values): 188 194 """ 189 195 Generate the distribution for parameter *name* given the parameter values … … 193 199 from the *pars* dictionary for parameter value and parameter dispersion. 194 200 """ 195 relative = name in model_info['partype']['pd-rel'] 196 limits = model_info['limits'][name] 197 disperser = pars.get(name+'_pd_type', 'gaussian') 198 value = pars.get(name, model_info['defaults'][name]) 199 npts = pars.get(name+'_pd_n', 0) 200 width = pars.get(name+'_pd', 0.0) 201 nsigma = pars.get(name+'_pd_nsigma', 3.0) 201 value = values.get(parameter.name, parameter.default) 202 relative = parameter.relative_pd 203 limits = parameter.limits 204 disperser = values.get(parameter.name+'_pd_type', 'gaussian') 205 npts = values.get(parameter.name+'_pd_n', 0) 206 width = values.get(parameter.name+'_pd', 0.0) 207 nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 208 if npts == 0 or width == 0: 209 return [value], [] 202 210 value, weight = weights.get_weights( 203 211 disperser, npts, width, nsigma, value, limits, relative) … … 220 228 def call_kernel(kernel, pars, cutoff=0, mono=False): 221 229 """ 222 Call *kernel* returned from :func:`make_kernel`with parameters *pars*.230 Call *kernel* returned from *model.make_kernel* with parameters *pars*. 223 231 224 232 *cutoff* is the limiting value for the product of dispersion weights used … … 228 236 with an error of about 1%, which is usually less than the measurement 229 237 uncertainty. 230 """ 231 fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 232 for name in kernel.fixed_pars] 238 239 *mono* is True if polydispersity should be set to none on all parameters. 240 """ 241 parameters = kernel.info['parameters'] 233 242 if mono: 234 pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 235 for name in kernel.pd_pars] 236 else: 237 pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 238 return kernel(fixed_pars, pd_pars, cutoff=cutoff) 243 active = lambda name: False 244 elif kernel.dim == '1d': 245 active = lambda name: name in parameters.pd_1d 246 elif kernel.dim == '2d': 247 active = lambda name: name in parameters.pd_2d 248 else: 249 active = lambda name: True 250 251 vw_pairs = [(get_weights(p, pars) if active(p.name) 252 else ([pars.get(p.name, p.default)], [])) 253 for p in parameters.call_parameters] 254 255 details, weights, values = build_details(kernel, vw_pairs) 256 return kernel(details, weights, values, cutoff) 257 258 def build_details(kernel, pairs): 259 values, weights = zip(*pairs) 260 if max([len(w) for w in weights]) > 1: 261 details = generate.poly_details(kernel.info, weights) 262 else: 263 details = kernel.info['mono_details'] 264 weights, values = [np.hstack(v) for v in (weights, values)] 265 weights = weights.astype(dtype=kernel.dtype) 266 values = values.astype(dtype=kernel.dtype) 267 return details, weights, values 268 239 269 240 270 def call_ER_VR(model_info, vol_pars): … … 259 289 260 290 261 def call_ER( info, pars):262 """ 263 Call the model ER function using * pars*.264 * info* is either *model.info* if you have a loaded model, or *kernel.info*265 if youhave a model kernel prepared for evaluation.266 """ 267 ER = info.get('ER', None)291 def call_ER(model_info, values): 292 """ 293 Call the model ER function using *values*. *model_info* is either 294 *model.info* if you have a loaded model, or *kernel.info* if you 295 have a model kernel prepared for evaluation. 296 """ 297 ER = model_info.get('ER', None) 268 298 if ER is None: 269 299 return 1.0 270 300 else: 271 vol_pars = [get_weights(info, pars, name) 272 for name in info['partype']['volume']] 301 vol_pars = [get_weights(parameter, values) 302 for parameter in model_info['parameters'] 303 if parameter.type == 'volume'] 273 304 value, weight = dispersion_mesh(vol_pars) 274 305 individual_radii = ER(*value) … … 276 307 return np.sum(weight*individual_radii) / np.sum(weight) 277 308 278 def call_VR( info, pars):309 def call_VR(model_info, values): 279 310 """ 280 311 Call the model VR function using *pars*. … … 282 313 if you have a model kernel prepared for evaluation. 283 314 """ 284 VR = info.get('VR', None)315 VR = model_info.get('VR', None) 285 316 if VR is None: 286 317 return 1.0 287 318 else: 288 vol_pars = [get_weights(info, pars, name) 289 for name in info['partype']['volume']] 319 vol_pars = [get_weights(parameter, values) 320 for parameter in model_info['parameters'] 321 if parameter.type == 'volume'] 290 322 value, weight = dispersion_mesh(vol_pars) 291 323 whole, part = VR(*value) -
sasmodels/direct_model.py
rea75043 rea75043 68 68 self.data_type = 'Iq' 69 69 70 partype = model.info['partype']71 72 70 if self.data_type == 'sesans': 73 71 … … 83 81 q_mono = sesans.make_all_q(data) 84 82 elif self.data_type == 'Iqxy': 85 if not partype['orientation'] and not partype['magnetic']:86 raise ValueError("not 2D without orientation or magnetic parameters")83 #if not model.info['parameters'].has_2d: 84 # raise ValueError("not 2D without orientation or magnetic parameters") 87 85 q = np.sqrt(data.qx_data**2 + data.qy_data**2) 88 86 qmin = getattr(data, 'qmin', 1e-16) -
sasmodels/generate.py
r9ef9dd9 rce896fd 21 21 22 22 *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 23 24 #define INVALID(v) (expr) returns False if v.parameter is invalid 25 for some parameter or other (e.g., v.bell_radius < v.radius). If 26 necessary, the expression can call a function. 23 27 24 28 These functions are defined in a kernel module .py script and an associated … … 65 69 The constructor code will generate the necessary vectors for computing 66 70 them with the desired polydispersity. 67 68 The available kernel parameters are defined as a list, with each parameter69 defined as a sublist with the following elements:70 71 *name* is the name that will be used in the call to the kernel72 function and the name that will be displayed to the user. Names73 should be lower case, with words separated by underscore. If74 acronyms are used, the whole acronym should be upper case.75 76 *units* should be one of *degrees* for angles, *Ang* for lengths,77 *1e-6/Ang^2* for SLDs.78 79 *default value* will be the initial value for the model when it80 is selected, or when an initial value is not otherwise specified.81 82 *limits = [lb, ub]* are the hard limits on the parameter value, used to83 limit the polydispersity density function. In the fit, the parameter limits84 given to the fit are the limits on the central value of the parameter.85 If there is polydispersity, it will evaluate parameter values outside86 the fit limits, but not outside the hard limits specified in the model.87 If there are no limits, use +/-inf imported from numpy.88 89 *type* indicates how the parameter will be used. "volume" parameters90 will be used in all functions. "orientation" parameters will be used91 in *Iqxy* and *Imagnetic*. "magnetic* parameters will be used in92 *Imagnetic* only. If *type* is the empty string, the parameter will93 be used in all of *Iq*, *Iqxy* and *Imagnetic*. "sld" parameters94 can automatically be promoted to magnetic parameters, each of which95 will have a magnitude and a direction, which may be different from96 other sld parameters.97 98 *description* is a short description of the parameter. This will99 be displayed in the parameter table and used as a tool tip for the100 parameter value in the user interface.101 102 71 The kernel module must set variables defining the kernel meta data: 103 72 … … 212 181 from __future__ import print_function 213 182 214 # TODO: identify model files which have changed since loading and reload them.215 216 import sys 183 #TODO: determine which functions are useful outside of generate 184 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 185 217 186 from os.path import abspath, dirname, join as joinpath, exists, basename, \ 218 splitext 187 splitext, getmtime 219 188 import re 220 189 import string 221 190 import warnings 222 from collections import namedtuple223 191 224 192 import numpy as np 225 193 226 PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 227 Parameter = namedtuple('Parameter', PARAMETER_FIELDS) 228 229 #TODO: determine which functions are useful outside of generate 230 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 231 232 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 194 from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo 195 196 # TODO: identify model files which have changed since loading and reload them. 197 198 TEMPLATE_ROOT = dirname(__file__) 233 199 234 200 F16 = np.dtype('float16') … … 239 205 except TypeError: 240 206 F128 = None 241 242 # Scale and background, which are parameters common to every form factor243 COMMON_PARAMETERS = [244 ["scale", "", 1, [0, np.inf], "", "Source intensity"],245 ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"],246 ]247 207 248 208 # Conversion from units defined in the parameter table for each model … … 338 298 raise ValueError("%r not found in %s" % (filename, search_path)) 339 299 300 340 301 def model_sources(model_info): 341 302 """ … … 346 307 return [_search(search_path, f) for f in model_info['source']] 347 308 348 # Pragmas for enable OpenCL features. Be sure to protect them so that they 349 # still compile even if OpenCL is not present. 350 _F16_PRAGMA = """\ 351 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 352 # pragma OPENCL EXTENSION cl_khr_fp16: enable 353 #endif 354 """ 355 356 _F64_PRAGMA = """\ 357 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 358 # pragma OPENCL EXTENSION cl_khr_fp64: enable 359 #endif 360 """ 309 def timestamp(model_info): 310 """ 311 Return a timestamp for the model corresponding to the most recently 312 changed file or dependency. 313 """ 314 source_files = (model_sources(model_info) 315 + model_templates() 316 + [model_info['filename']]) 317 newest = max(getmtime(f) for f in source_files) 318 return newest 361 319 362 320 def convert_type(source, dtype): … … 369 327 if dtype == F16: 370 328 fbytes = 2 371 source = _ F16_PRAGMA + _convert_type(source, "half", "f")329 source = _convert_type(source, "float", "f") 372 330 elif dtype == F32: 373 331 fbytes = 4 … … 375 333 elif dtype == F64: 376 334 fbytes = 8 377 source = _F64_PRAGMA + source # Sourceis already double335 # no need to convert if it is already double 378 336 elif dtype == F128: 379 337 fbytes = 16 … … 418 376 419 377 420 LOOP_OPEN = """\ 421 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 422 const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 423 const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 378 _template_cache = {} 379 def load_template(filename): 380 path = joinpath(TEMPLATE_ROOT, filename) 381 mtime = getmtime(path) 382 if filename not in _template_cache or mtime > _template_cache[filename][0]: 383 with open(path) as fid: 384 _template_cache[filename] = (mtime, fid.read(), path) 385 return _template_cache[filename][1] 386 387 def model_templates(): 388 # TODO: fails DRY; templates are listed in two places. 389 # should instead have model_info contain a list of paths 390 return [joinpath(TEMPLATE_ROOT, filename) 391 for filename in ('kernel_header.c', 'kernel_iq.c')] 392 393 394 _FN_TEMPLATE = """\ 395 double %(name)s(%(pars)s); 396 double %(name)s(%(pars)s) { 397 %(body)s 398 } 399 400 424 401 """ 425 def build_polydispersity_loops(pd_pars): 426 """ 427 Build polydispersity loops 428 429 Returns loop opening and loop closing 430 """ 431 depth = 4 432 offset = "" 433 loop_head = [] 434 loop_end = [] 435 for name in pd_pars: 436 subst = {'name': name, 'offset': offset} 437 loop_head.append(indent(LOOP_OPEN % subst, depth)) 438 loop_end.insert(0, (" "*depth) + "}") 439 offset += '+N' + name 440 depth += 2 441 return "\n".join(loop_head), "\n".join(loop_end) 442 443 C_KERNEL_TEMPLATE = None 402 403 def _gen_fn(name, pars, body): 404 """ 405 Generate a function given pars and body. 406 407 Returns the following string:: 408 409 double fn(double a, double b, ...); 410 double fn(double a, double b, ...) { 411 .... 412 } 413 """ 414 par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 415 return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 416 417 def _call_pars(prefix, pars): 418 """ 419 Return a list of *prefix.parameter* from parameter items. 420 """ 421 return [p.as_call_reference(prefix) for p in pars] 422 423 _IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 424 flags=re.MULTILINE) 425 def _have_Iqxy(sources): 426 """ 427 Return true if any file defines Iqxy. 428 429 Note this is not a C parser, and so can be easily confused by 430 non-standard syntax. Also, it will incorrectly identify the following 431 as having Iqxy:: 432 433 /* 434 double Iqxy(qx, qy, ...) { ... fill this in later ... } 435 */ 436 437 If you want to comment out an Iqxy function, use // on the front of the 438 line instead. 439 """ 440 for code in sources: 441 if _IQXY_PATTERN.search(code): 442 return True 443 else: 444 return False 445 444 446 def make_source(model_info): 445 447 """ … … 460 462 # for computing volume even if we allow non-disperse volume parameters. 461 463 462 # Load template 463 global C_KERNEL_TEMPLATE 464 if C_KERNEL_TEMPLATE is None: 465 with open(C_KERNEL_TEMPLATE_PATH) as fid: 466 C_KERNEL_TEMPLATE = fid.read() 467 468 # Load additional sources 469 source = [open(f).read() for f in model_sources(model_info)] 470 471 # Prepare defines 472 defines = [] 473 partype = model_info['partype'] 474 pd_1d = partype['pd-1d'] 475 pd_2d = partype['pd-2d'] 476 fixed_1d = partype['fixed-1d'] 477 fixed_2d = partype['fixed-1d'] 478 479 iq_parameters = [p.name 480 for p in model_info['parameters'][2:] # skip scale, background 481 if p.name in set(fixed_1d + pd_1d)] 482 iqxy_parameters = [p.name 483 for p in model_info['parameters'][2:] # skip scale, background 484 if p.name in set(fixed_2d + pd_2d)] 485 volume_parameters = [p.name 486 for p in model_info['parameters'] 487 if p.type == 'volume'] 488 489 # Fill in defintions for volume parameters 490 if volume_parameters: 491 defines.append(('VOLUME_PARAMETERS', 492 ','.join(volume_parameters))) 493 defines.append(('VOLUME_WEIGHT_PRODUCT', 494 '*'.join(p + '_w' for p in volume_parameters))) 495 496 # Generate form_volume function from body only 464 partable = model_info['parameters'] 465 466 # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 467 # Note that scale and volume are not possible types. 468 469 # Load templates and user code 470 kernel_header = load_template('kernel_header.c') 471 kernel_code = load_template('kernel_iq.c') 472 user_code = [open(f).read() for f in model_sources(model_info)] 473 474 # Build initial sources 475 source = [kernel_header] + user_code 476 477 # Make parameters for q, qx, qy so that we can use them in declarations 478 q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 479 # Generate form_volume function, etc. from body only 497 480 if model_info['form_volume'] is not None: 498 if volume_parameters: 499 vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 500 else: 501 vol_par_decl = 'void' 502 defines.append(('VOLUME_PARAMETER_DECLARATIONS', 503 vol_par_decl)) 504 fn = """\ 505 double form_volume(VOLUME_PARAMETER_DECLARATIONS); 506 double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 507 %(body)s 508 } 509 """ % {'body':model_info['form_volume']} 510 source.append(fn) 511 512 # Fill in definitions for Iq parameters 513 defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq')) 514 defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 515 if fixed_1d: 516 defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 517 ', \\\n '.join('const double %s' % p for p in fixed_1d))) 518 if pd_1d: 519 defines.append(('IQ_WEIGHT_PRODUCT', 520 '*'.join(p + '_w' for p in pd_1d))) 521 defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 522 ', \\\n '.join('const int N%s' % p for p in pd_1d))) 523 defines.append(('IQ_DISPERSION_LENGTH_SUM', 524 '+'.join('N' + p for p in pd_1d))) 525 open_loops, close_loops = build_polydispersity_loops(pd_1d) 526 defines.append(('IQ_OPEN_LOOPS', 527 open_loops.replace('\n', ' \\\n'))) 528 defines.append(('IQ_CLOSE_LOOPS', 529 close_loops.replace('\n', ' \\\n'))) 481 pars = partable.form_volume_parameters 482 source.append(_gen_fn('form_volume', pars, model_info['form_volume'])) 530 483 if model_info['Iq'] is not None: 531 defines.append(('IQ_PARAMETER_DECLARATIONS', 532 ', '.join('double ' + p for p in iq_parameters))) 533 fn = """\ 534 double Iq(double q, IQ_PARAMETER_DECLARATIONS); 535 double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 536 %(body)s 537 } 538 """ % {'body':model_info['Iq']} 539 source.append(fn) 540 541 # Fill in definitions for Iqxy parameters 542 defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 543 defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 544 if fixed_2d: 545 defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 546 ', \\\n '.join('const double %s' % p for p in fixed_2d))) 547 if pd_2d: 548 defines.append(('IQXY_WEIGHT_PRODUCT', 549 '*'.join(p + '_w' for p in pd_2d))) 550 defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 551 ', \\\n '.join('const int N%s' % p for p in pd_2d))) 552 defines.append(('IQXY_DISPERSION_LENGTH_SUM', 553 '+'.join('N' + p for p in pd_2d))) 554 open_loops, close_loops = build_polydispersity_loops(pd_2d) 555 defines.append(('IQXY_OPEN_LOOPS', 556 open_loops.replace('\n', ' \\\n'))) 557 defines.append(('IQXY_CLOSE_LOOPS', 558 close_loops.replace('\n', ' \\\n'))) 484 pars = [q] + partable.iq_parameters 485 source.append(_gen_fn('Iq', pars, model_info['Iq'])) 559 486 if model_info['Iqxy'] is not None: 560 defines.append(('IQXY_PARAMETER_DECLARATIONS', 561 ', '.join('double ' + p for p in iqxy_parameters))) 562 fn = """\ 563 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 564 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 565 %(body)s 566 } 567 """ % {'body':model_info['Iqxy']} 568 source.append(fn) 569 570 # Need to know if we have a theta parameter for Iqxy; it is not there 571 # for the magnetic sphere model, for example, which has a magnetic 572 # orientation but no shape orientation. 573 if 'theta' in pd_2d: 574 defines.append(('IQXY_HAS_THETA', '1')) 575 576 #for d in defines: print(d) 577 defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 578 sources = '\n\n'.join(source) 579 return C_KERNEL_TEMPLATE % { 580 'DEFINES': defines, 581 'SOURCES': sources, 582 } 583 584 def categorize_parameters(pars): 585 """ 586 Build parameter categories out of the the parameter definitions. 587 588 Returns a dictionary of categories. 589 590 Note: these categories are subject to change, depending on the needs of 591 the UI and the needs of the kernel calling function. 592 593 The categories are as follows: 594 595 * *volume* list of volume parameter names 596 * *orientation* list of orientation parameters 597 * *magnetic* list of magnetic parameters 598 * *<empty string>* list of parameters that have no type info 599 600 Each parameter is in one and only one category. 601 602 The following derived categories are created: 603 604 * *fixed-1d* list of non-polydisperse parameters for 1D models 605 * *pd-1d* list of polydisperse parameters for 1D models 606 * *fixed-2d* list of non-polydisperse parameters for 2D models 607 * *pd-d2* list of polydisperse parameters for 2D models 608 """ 609 partype = { 610 'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 611 'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 612 'pd-rel': set(), 613 } 614 615 for p in pars: 616 if p.type == 'volume': 617 partype['pd-1d'].append(p.name) 618 partype['pd-2d'].append(p.name) 619 partype['pd-rel'].add(p.name) 620 elif p.type == 'magnetic': 621 partype['fixed-2d'].append(p.name) 622 elif p.type == 'orientation': 623 partype['pd-2d'].append(p.name) 624 elif p.type in ('', 'sld'): 625 partype['fixed-1d'].append(p.name) 626 partype['fixed-2d'].append(p.name) 627 else: 628 raise ValueError("unknown parameter type %r" % p.type) 629 partype[p.type].append(p.name) 630 631 return partype 632 633 def process_parameters(model_info): 634 """ 635 Process parameter block, precalculating parameter details. 636 """ 637 # convert parameters into named tuples 638 for p in model_info['parameters']: 639 if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 640 p[4] = 'sld' 641 # TODO: make sure all models explicitly label their sld parameters 642 #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 643 644 pars = [Parameter(*p) for p in model_info['parameters']] 645 # Fill in the derived attributes 646 model_info['parameters'] = pars 647 partype = categorize_parameters(pars) 648 model_info['limits'] = dict((p.name, p.limits) for p in pars) 649 model_info['partype'] = partype 650 model_info['defaults'] = dict((p.name, p.default) for p in pars) 651 if model_info.get('demo', None) is None: 652 model_info['demo'] = model_info['defaults'] 653 model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 487 pars = [qx, qy] + partable.iqxy_parameters 488 source.append(_gen_fn('Iqxy', pars, model_info['Iqxy'])) 489 490 # Define the parameter table 491 source.append("#define PARAMETER_TABLE \\") 492 source.append("\\\n".join(p.as_definition() 493 for p in partable.kernel_parameters)) 494 495 # Define the function calls 496 if partable.form_volume_parameters: 497 refs = _call_pars("v.", partable.form_volume_parameters) 498 call_volume = "#define CALL_VOLUME(v) form_volume(%s)" % (",".join(refs)) 499 else: 500 # Model doesn't have volume. We could make the kernel run a little 501 # faster by not using/transferring the volume normalizations, but 502 # the ifdef's reduce readability more than is worthwhile. 503 call_volume = "#define CALL_VOLUME(v) 1.0" 504 source.append(call_volume) 505 506 refs = ["q[i]"] + _call_pars("v.", partable.iq_parameters) 507 call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs)) 508 if _have_Iqxy(user_code): 509 # Call 2D model 510 refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v.", partable.iqxy_parameters) 511 call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs)) 512 else: 513 # Call 1D model with sqrt(qx^2 + qy^2) 514 warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 515 # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 516 pars_sqrt = ["sqrt(q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:] 517 call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt)) 518 519 # Fill in definitions for numbers of parameters 520 source.append("#define MAX_PD %s"%partable.max_pd) 521 source.append("#define NPARS %d"%partable.npars) 522 523 # TODO: allow mixed python/opencl kernels? 524 525 # define the Iq kernel 526 source.append("#define KERNEL_NAME %s_Iq"%model_info['name']) 527 source.append(call_iq) 528 source.append(kernel_code) 529 source.append("#undef CALL_IQ") 530 source.append("#undef KERNEL_NAME") 531 532 # define the Iqxy kernel from the same source with different #defines 533 source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name']) 534 source.append(call_iqxy) 535 source.append(kernel_code) 536 source.append("#undef CALL_IQ") 537 source.append("#undef KERNEL_NAME") 538 539 return '\n'.join(source) 540 541 class CoordinationDetails(object): 542 def __init__(self, model_info): 543 parameters = model_info['parameters'] 544 max_pd = parameters.max_pd 545 npars = parameters.npars 546 par_offset = 4*max_pd 547 self.details = np.zeros(par_offset + 3*npars + 4, 'i4') 548 549 # generate views on different parts of the array 550 self._pd_par = self.details[0*max_pd:1*max_pd] 551 self._pd_length = self.details[1*max_pd:2*max_pd] 552 self._pd_offset = self.details[2*max_pd:3*max_pd] 553 self._pd_stride = self.details[3*max_pd:4*max_pd] 554 self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars] 555 self._par_coord = self.details[par_offset+1*npars:par_offset+2*npars] 556 self._pd_coord = self.details[par_offset+2*npars:par_offset+3*npars] 557 558 # theta_par is fixed 559 self.details[-1] = parameters.theta_offset 560 561 @property 562 def ctypes(self): return self.details.ctypes 563 @property 564 def pd_par(self): return self._pd_par 565 @property 566 def pd_length(self): return self._pd_length 567 @property 568 def pd_offset(self): return self._pd_offset 569 @property 570 def pd_stride(self): return self._pd_stride 571 @property 572 def pd_coord(self): return self._pd_coord 573 @property 574 def par_coord(self): return self._par_coord 575 @property 576 def par_offset(self): return self._par_offset 577 @property 578 def num_coord(self): return self.details[-2] 579 @num_coord.setter 580 def num_coord(self, v): self.details[-2] = v 581 @property 582 def total_pd(self): return self.details[-3] 583 @total_pd.setter 584 def total_pd(self, v): self.details[-3] = v 585 @property 586 def num_active(self): return self.details[-4] 587 @num_active.setter 588 def num_active(self, v): self.details[-4] = v 589 590 def show(self): 591 print("total_pd", self.total_pd) 592 print("num_active", self.num_active) 593 print("pd_par", self.pd_par) 594 print("pd_length", self.pd_length) 595 print("pd_offset", self.pd_offset) 596 print("pd_stride", self.pd_stride) 597 print("par_offsets", self.par_offset) 598 print("num_coord", self.num_coord) 599 print("par_coord", self.par_coord) 600 print("pd_coord", self.pd_coord) 601 print("theta par", self.details[-1]) 602 603 def mono_details(model_info): 604 details = CoordinationDetails(model_info) 605 # The zero defaults for monodisperse systems are mostly fine 606 details.par_offset[:] = np.arange(2, len(details.par_offset)+2) 607 return details 608 609 def poly_details(model_info, weights): 610 #print("weights",weights) 611 weights = weights[2:] # Skip scale and background 612 613 # Decreasing list of polydispersity lengths 614 # Note: the reversing view, x[::-1], does not require a copy 615 pd_length = np.array([len(w) for w in weights]) 616 num_active = np.sum(pd_length>1) 617 if num_active > model_info['parameters'].max_pd: 618 raise ValueError("Too many polydisperse parameters") 619 620 pd_offset = np.cumsum(np.hstack((0, pd_length))) 621 idx = np.argsort(pd_length)[::-1][:num_active] 622 par_length = np.array([max(len(w),1) for w in weights]) 623 pd_stride = np.cumprod(np.hstack((1, par_length[idx]))) 624 par_offsets = np.cumsum(np.hstack((2, par_length))) 625 626 details = CoordinationDetails(model_info) 627 details.pd_par[:num_active] = idx 628 details.pd_length[:num_active] = pd_length[idx] 629 details.pd_offset[:num_active] = pd_offset[idx] 630 details.pd_stride[:num_active] = pd_stride[:-1] 631 details.par_offset[:] = par_offsets[:-1] 632 details.total_pd = pd_stride[-1] 633 details.num_active = num_active 634 # Without constraints coordinated parameters are just the pd parameters 635 details.par_coord[:num_active] = idx 636 details.pd_coord[:num_active] = 2**np.arange(num_active) 637 details.num_coord = num_active 638 #details.show() 639 return details 640 641 def constrained_poly_details(model_info, weights, constraints): 642 # Need to find the independently varying pars and sort them 643 # Need to build a coordination list for the dependent variables 644 # Need to generate a constraints function which takes values 645 # and weights, returning par blocks 646 raise NotImplementedError("Can't handle constraints yet") 647 648 649 def create_default_functions(model_info): 650 """ 651 Autogenerate missing functions, such as Iqxy from Iq. 652 653 This only works for Iqxy when Iq is written in python. :func:`make_source` 654 performs a similar role for Iq written in C. 655 """ 656 if callable(model_info['Iq']) and model_info['Iqxy'] is None: 657 partable = model_info['parameters'] 658 if partable.has_2d: 659 raise ValueError("Iqxy model is missing") 660 Iq = model_info['Iq'] 661 def Iqxy(qx, qy, **kw): 662 return Iq(np.sqrt(qx**2 + qy**2), **kw) 663 model_info['Iqxy'] = Iqxy 664 654 665 655 666 def make_model_info(kernel_module): … … 678 689 model variants (e.g., the list of cases) or is None if there are no 679 690 model variants 680 * *defaults* is the *{parameter: value}* table built from the parameter 681 description table. 682 * *limits* is the *{parameter: [min, max]}* table built from the 683 parameter description table. 684 * *partypes* categorizes the model parameters. See 691 * *par_type* categorizes the model parameters. See 685 692 :func:`categorize_parameters` for details. 686 693 * *demo* contains the *{parameter: value}* map used in compare (and maybe … … 699 706 *model_info* blocks for the composition objects. This allows us to 700 707 build complete product and mixture models from just the info. 701 702 708 """ 703 709 # TODO: maybe turn model_info into a class ModelDefinition 704 parameters = COMMON_PARAMETERS + kernel_module.parameters 710 #print("make parameter table", kernel_module.parameters) 711 parameters = make_parameter_table(kernel_module.parameters) 705 712 filename = abspath(kernel_module.__file__) 706 713 kernel_id = splitext(basename(filename))[0] … … 712 719 filename=abspath(kernel_module.__file__), 713 720 name=name, 714 title= kernel_module.title,715 description= kernel_module.description,721 title=getattr(kernel_module, 'title', name+" model"), 722 description=getattr(kernel_module, 'description', 'no description'), 716 723 parameters=parameters, 717 724 composition=None, … … 720 727 single=getattr(kernel_module, 'single', True), 721 728 structure_factor=getattr(kernel_module, 'structure_factor', False), 729 profile_axes=getattr(kernel_module, 'profile_axes', ['x','y']), 722 730 variant_info=getattr(kernel_module, 'invariant_info', None), 723 731 demo=getattr(kernel_module, 'demo', None), … … 727 735 tests=getattr(kernel_module, 'tests', []), 728 736 ) 729 process_parameters(model_info)737 set_demo(model_info, getattr(kernel_module, 'demo', None)) 730 738 # Check for optional functions 731 functions = "ER VR form_volume Iq Iqxy shape sesans".split()739 functions = "ER VR form_volume Iq Iqxy profile sesans".split() 732 740 model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 741 create_default_functions(model_info) 742 # Precalculate the monodisperse parameters 743 # TODO: make this a lazy evaluator 744 # make_model_info is called for every model on sasview startup 745 model_info['mono_details'] = mono_details(model_info) 733 746 return model_info 734 747 … … 782 795 783 796 784 785 797 def demo_time(): 786 798 """ … … 798 810 Program which prints the source produced by the model. 799 811 """ 812 import sys 813 from sasmodels.core import make_model_by_name 800 814 if len(sys.argv) <= 1: 801 815 print("usage: python -m sasmodels.generate modelname") 802 816 else: 803 817 name = sys.argv[1] 804 import sasmodels.models 805 __import__('sasmodels.models.' + name) 806 model = getattr(sasmodels.models, name) 807 model_info = make_model_info(model) 818 model_info = make_model_by_name(name) 808 819 source = make_source(model_info) 809 820 print(source) -
sasmodels/kernelcl.py
r8e0d974 rba32cdd 48 48 harmless, albeit annoying. 49 49 """ 50 from __future__ import print_function 50 51 import os 51 52 import warnings … … 73 74 # of polydisperse parameters. 74 75 MAX_LOOPS = 2048 76 77 78 # Pragmas for enable OpenCL features. Be sure to protect them so that they 79 # still compile even if OpenCL is not present. 80 _F16_PRAGMA = """\ 81 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 82 # pragma OPENCL EXTENSION cl_khr_fp16: enable 83 #endif 84 """ 85 86 _F64_PRAGMA = """\ 87 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 88 # pragma OPENCL EXTENSION cl_khr_fp64: enable 89 #endif 90 """ 75 91 76 92 … … 142 158 raise RuntimeError("%s not supported for devices"%dtype) 143 159 144 source = generate.convert_type(source, dtype) 160 source_list = [generate.convert_type(source, dtype)] 161 162 if dtype == generate.F16: 163 source_list.insert(0, _F16_PRAGMA) 164 elif dtype == generate.F64: 165 source_list.insert(0, _F64_PRAGMA) 166 145 167 # Note: USE_SINCOS makes the intel cpu slower under opencl 146 168 if context.devices[0].type == cl.device_type.GPU: 147 source = "#define USE_SINCOS\n" + source169 source_list.insert(0, "#define USE_SINCOS\n") 148 170 options = (get_fast_inaccurate_build_options(context.devices[0]) 149 171 if fast else []) 172 source = "\n".join(source_list) 150 173 program = cl.Program(context, source).build(options=options) 151 174 return program … … 220 243 key = "%s-%s-%s"%(name, dtype, fast) 221 244 if key not in self.compiled: 222 #print("compiling",name)245 print("compiling",name) 223 246 dtype = np.dtype(dtype) 224 program = compile_model(self.get_context(dtype), source, dtype, fast) 247 program = compile_model(self.get_context(dtype), 248 str(source), dtype, fast) 225 249 self.compiled[key] = program 226 250 return self.compiled[key] … … 314 338 self.program = None 315 339 316 def __call__(self, q_vectors):340 def make_kernel(self, q_vectors): 317 341 if self.program is None: 318 342 compiler = environment().compile_program 319 self.program = compiler(self.info['name'], self.source, self.dtype,320 self. fast)343 self.program = compiler(self.info['name'], self.source, 344 self.dtype, self.fast) 321 345 is_2d = len(q_vectors) == 2 322 346 kernel_name = generate.kernel_name(self.info, is_2d) … … 365 389 # at this point, so instead using 32, which is good on the set of 366 390 # architectures tested so far. 367 self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 391 if self.is_2d: 392 # Note: 17 rather than 15 because results is 2 elements 393 # longer than input. 394 width = ((self.nq+17)//16)*16 395 self.q = np.empty((width, 2), dtype=dtype) 396 self.q[:self.nq, 0] = q_vectors[0] 397 self.q[:self.nq, 1] = q_vectors[1] 398 else: 399 # Note: 33 rather than 31 because results is 2 elements 400 # longer than input. 401 width = ((self.nq+33)//32)*32 402 self.q = np.empty(width, dtype=dtype) 403 self.q[:self.nq] = q_vectors[0] 404 self.global_size = [self.q.shape[0]] 368 405 context = env.get_context(self.dtype) 369 self.global_size = [self.q_vectors[0].size]370 406 #print("creating inputs of size", self.global_size) 371 self.q_buffers = [ 372 cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 373 for q in self.q_vectors 374 ] 407 self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 408 hostbuf=self.q) 375 409 376 410 def release(self): … … 378 412 Free the memory. 379 413 """ 380 for b in self.q_buffers:381 b.release()382 self.q_buffers = []414 if self.q is not None: 415 self.q.release() 416 self.q = None 383 417 384 418 def __del__(self): … … 406 440 """ 407 441 def __init__(self, kernel, model_info, q_vectors, dtype): 442 max_pd = model_info['max_pd'] 443 npars = len(model_info['parameters'])-2 408 444 q_input = GpuInput(q_vectors, dtype) 445 self.dtype = dtype 446 self.dim = '2d' if q_input.is_2d else '1d' 409 447 self.kernel = kernel 410 448 self.info = model_info 411 self.res = np.empty(q_input.nq, q_input.dtype) 412 dim = '2d' if q_input.is_2d else '1d' 413 self.fixed_pars = model_info['partype']['fixed-' + dim] 414 self.pd_pars = model_info['partype']['pd-' + dim] 449 self.pd_stop_index = 4*max_pd-1 450 # plus three for the normalization values 451 self.result = np.empty(q_input.nq+3, q_input.dtype) 415 452 416 453 # Inputs and outputs for each kernel call … … 418 455 env = environment() 419 456 self.queue = env.get_queue(dtype) 420 self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 421 2 * MAX_LOOPS * q_input.dtype.itemsize) 422 self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 457 458 # details is int32 data, padded to an 8 integer boundary 459 size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 460 self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 423 461 q_input.global_size[0] * q_input.dtype.itemsize) 424 self.q_input = q_input 425 426 self._need_release = [ self.loops_b, self.res_b, self.q_input]427 428 def __call__(self, fixed_pars, pd_pars, cutoff):462 self.q_input = q_input # allocated by GpuInput above 463 464 self._need_release = [ self.result_b, self.q_input ] 465 466 def __call__(self, details, weights, values, cutoff): 429 467 real = (np.float32 if self.q_input.dtype == generate.F32 430 468 else np.float64 if self.q_input.dtype == generate.F64 431 469 else np.float16 if self.q_input.dtype == generate.F16 432 470 else np.float32) # will never get here, so use np.float32 433 434 #print "pars", fixed_pars, pd_pars 435 res_bi = self.res_b 436 nq = np.uint32(self.q_input.nq) 437 if pd_pars: 438 cutoff = real(cutoff) 439 loops_N = [np.uint32(len(p[0])) for p in pd_pars] 440 loops = np.hstack(pd_pars) \ 441 if pd_pars else np.empty(0, dtype=self.q_input.dtype) 442 loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 443 #print("loops",Nloops, loops) 444 445 #import sys; print("opencl eval",pars) 446 #print("opencl eval",pars) 447 if len(loops) > 2 * MAX_LOOPS: 448 raise ValueError("too many polydispersity points") 449 450 loops_bi = self.loops_b 451 cl.enqueue_copy(self.queue, loops_bi, loops) 452 loops_l = cl.LocalMemory(len(loops.data)) 453 #ctx = environment().context 454 #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops) 455 dispersed = [loops_bi, loops_l, cutoff] + loops_N 456 else: 457 dispersed = [] 458 fixed = [real(p) for p in fixed_pars] 459 args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed 471 assert details.dtype == np.int32 472 assert weights.dtype == real and values.dtype == real 473 474 context = self.queue.context 475 details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 476 hostbuf=details) 477 weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 478 hostbuf=weights) 479 values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 480 hostbuf=values) 481 482 start, stop = 0, self.details[self.pd_stop_index] 483 args = [ 484 np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 485 self.details_b, self.weights_b, self.values_b, 486 self.q_input.q_b, self.result_b, real(cutoff), 487 ] 460 488 self.kernel(self.queue, self.q_input.global_size, None, *args) 461 cl.enqueue_copy(self.queue, self.res, res_bi) 462 463 return self.res 489 cl.enqueue_copy(self.queue, self.result, self.result_b) 490 [v.release() for v in details_b, weights_b, values_b] 491 492 return self.result[:self.nq] 464 493 465 494 def release(self): -
sasmodels/kerneldll.py
r6ad0e87 rd19962c 50 50 import tempfile 51 51 import ctypes as ct 52 from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float 53 import _ctypes 52 from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float 54 53 55 54 import numpy as np … … 81 80 COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 82 81 if "SAS_OPENMP" in os.environ: 83 COMPILE = COMPILE +" -fopenmp"82 COMPILE += " -fopenmp" 84 83 else: 85 84 COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" … … 90 89 91 90 92 def dll_ path(model_info, dtype="double"):93 """ 94 Path to the compiled model defined by *model_info*.95 """96 from os.path import join as joinpath, split as splitpath, splitext97 b asename = splitext(splitpath(model_info['filename'])[1])[0]98 if np.dtype(dtype) == generate.F32:99 basename += "32" 100 elif np.dtype(dtype) == generate.F64:101 basename += "64"102 else:103 basename += "128"104 return joinpath(DLL_PATH, basename+'.so')105 91 def dll_name(model_info, dtype): 92 """ 93 Name of the dll containing the model. This is the base file name without 94 any path or extension, with a form such as 'sas_sphere32'. 95 """ 96 bits = 8*dtype.itemsize 97 return "sas_%s%d"%(model_info['id'], bits) 98 99 def dll_path(model_info, dtype): 100 """ 101 Complete path to the dll for the model. Note that the dll may not 102 exist yet if it hasn't been compiled. 103 """ 104 return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so") 106 105 107 106 def make_dll(source, model_info, dtype="double"): … … 133 132 dtype = generate.F64 # Force 64-bit dll 134 133 135 if dtype == generate.F32: # 32-bit dll136 tempfile_prefix = 'sas_' + model_info['name'] + '32_'137 elif dtype == generate.F64:138 tempfile_prefix = 'sas_' + model_info['name'] + '64_'139 else:140 tempfile_prefix = 'sas_' + model_info['name'] + '128_'141 142 134 source = generate.convert_type(source, dtype) 143 source_files = generate.model_sources(model_info) + [model_info['filename']]135 newest = generate.timestamp(model_info) 144 136 dll = dll_path(model_info, dtype) 145 newest = max(os.path.getmtime(f) for f in source_files)146 137 if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 147 # Replace with a proper temp file148 fid, filename = tempfile.mkstemp(suffix=".c", prefix= tempfile_prefix)138 basename = dll_name(model_info, dtype) + "_" 139 fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 149 140 os.fdopen(fid, "w").write(source) 150 141 command = COMPILE%{"source":filename, "output":dll} … … 171 162 return DllModel(filename, model_info, dtype=dtype) 172 163 173 174 IQ_ARGS = [c_void_p, c_void_p, c_int]175 IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int]176 177 164 class DllModel(object): 178 165 """ … … 197 184 198 185 def _load_dll(self): 199 Nfixed1d = len(self.info['partype']['fixed-1d'])200 Nfixed2d = len(self.info['partype']['fixed-2d'])201 Npd1d = len(self.info['partype']['pd-1d'])202 Npd2d = len(self.info['partype']['pd-2d'])203 204 186 #print("dll", self.dllpath) 205 187 try: … … 212 194 else c_double if self.dtype == generate.F64 213 195 else c_longdouble) 214 pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 215 pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 196 197 # int, int, int, int*, double*, double*, double*, double*, double*, double 198 argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 216 199 self.Iq = self.dll[generate.kernel_name(self.info, False)] 217 self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d218 219 200 self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 220 self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d 221 222 self.release() 201 self.Iq.argtypes = argtypes 202 self.Iqxy.argtypes = argtypes 223 203 224 204 def __getstate__(self): … … 229 209 self.dll = None 230 210 231 def __call__(self, q_vectors):211 def make_kernel(self, q_vectors): 232 212 q_input = PyInput(q_vectors, self.dtype) 233 213 if self.dll is None: self._load_dll() … … 272 252 """ 273 253 def __init__(self, kernel, model_info, q_input): 254 self.kernel = kernel 274 255 self.info = model_info 275 256 self.q_input = q_input 276 self.kernel = kernel 277 self.res = np.empty(q_input.nq, q_input.dtype) 278 dim = '2d' if q_input.is_2d else '1d' 279 self.fixed_pars = model_info['partype']['fixed-' + dim] 280 self.pd_pars = model_info['partype']['pd-' + dim] 281 282 # In dll kernel, but not in opencl kernel 283 self.p_res = self.res.ctypes.data 284 285 def __call__(self, fixed_pars, pd_pars, cutoff): 257 self.dtype = q_input.dtype 258 self.dim = '2d' if q_input.is_2d else '1d' 259 self.result = np.empty(q_input.nq+3, q_input.dtype) 260 261 def __call__(self, details, weights, values, cutoff): 286 262 real = (np.float32 if self.q_input.dtype == generate.F32 287 263 else np.float64 if self.q_input.dtype == generate.F64 288 264 else np.float128) 289 290 nq = c_int(self.q_input.nq) 291 if pd_pars: 292 cutoff = real(cutoff) 293 loops_N = [np.uint32(len(p[0])) for p in pd_pars] 294 loops = np.hstack(pd_pars) 295 loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 296 p_loops = loops.ctypes.data 297 dispersed = [p_loops, cutoff] + loops_N 298 else: 299 dispersed = [] 300 fixed = [real(p) for p in fixed_pars] 301 args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed 302 #print(pars) 265 assert isinstance(details, generate.CoordinationDetails) 266 assert weights.dtype == real and values.dtype == real 267 268 start, stop = 0, details.total_pd 269 #print("in kerneldll") 270 #print("weights", weights) 271 #print("values", values) 272 args = [ 273 self.q_input.nq, # nq 274 start, # pd_start 275 stop, # pd_stop pd_stride[MAX_PD] 276 details.ctypes.data, # problem 277 weights.ctypes.data, # weights 278 values.ctypes.data, #pars 279 self.q_input.q.ctypes.data, #q 280 self.result.ctypes.data, # results 281 real(cutoff), # cutoff 282 ] 303 283 self.kernel(*args) 304 305 return self.res 284 return self.result[:-3] 306 285 307 286 def release(self): -
sasmodels/kernelpy.py
ra84a0ca r48fbd50 53 53 self.dtype = dtype 54 54 self.is_2d = (len(q_vectors) == 2) 55 self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 56 self.q_pointers = [q.ctypes.data for q in self.q_vectors] 55 if self.is_2d: 56 self.q = np.empty((self.nq, 2), dtype=dtype) 57 self.q[:, 0] = q_vectors[0] 58 self.q[:, 1] = q_vectors[1] 59 else: 60 self.q = np.empty(self.nq, dtype=dtype) 61 self.q[:self.nq] = q_vectors[0] 57 62 58 63 def release(self): … … 60 65 Free resources associated with the model inputs. 61 66 """ 62 self.q _vectors = []67 self.q = None 63 68 64 69 class PyKernel(object): -
sasmodels/mixture.py
r72a081d r69aa451 14 14 import numpy as np 15 15 16 from . generate import process_parameters, COMMON_PARAMETERS, Parameter16 from .modelinfo import Parameter, ParameterTable 17 17 18 18 SCALE=0 … … 34 34 35 35 # Build new parameter list 36 pars = COMMON_PARAMETERS +[]36 pars = [] 37 37 for k, part in enumerate(parts): 38 38 # Parameter prefix per model, A_, B_, ... 39 # Note that prefix must also be applied to id and length_control 40 # to support vector parameters 39 41 prefix = chr(ord('A')+k) + '_' 40 for p in part['parameters']: 41 # No background on the individual mixture elements 42 if p.name == 'background': 43 continue 44 # TODO: promote Parameter to a full class 45 # this code knows too much about the implementation! 46 p = list(p) 47 if p[0] == 'scale': # allow model subtraction 48 p[3] = [-np.inf, np.inf] # limits 49 p[0] = prefix+p[0] # relabel parameters with A_, ... 42 pars.append(Parameter(prefix+'scale')) 43 for p in part['parameters'].kernel_pars: 44 p = copy(p) 45 p.name = prefix+p.name 46 p.id = prefix+p.id 47 if p.length_control is not None: 48 p.length_control = prefix+p.length_control 50 49 pars.append(p) 50 partable = ParameterTable(pars) 51 51 52 52 model_info = {} … … 58 58 model_info['docs'] = model_info['title'] 59 59 model_info['category'] = "custom" 60 model_info['parameters'] = par s60 model_info['parameters'] = partable 61 61 #model_info['single'] = any(part['single'] for part in parts) 62 62 model_info['structure_factor'] = False … … 67 67 # Remember the component info blocks so we can build the model 68 68 model_info['composition'] = ('mixture', parts) 69 process_parameters(model_info)70 return model_info71 69 72 70 -
sasmodels/model_test.py
ra84a0ca r48fbd50 51 51 52 52 from .core import list_models, load_model_info, build_model, HAVE_OPENCL 53 from .core import make_kernel,call_kernel, call_ER, call_VR53 from .core import call_kernel, call_ER, call_VR 54 54 from .exception import annotate_exception 55 55 … … 187 187 Qx, Qy = zip(*x) 188 188 q_vectors = [np.array(Qx), np.array(Qy)] 189 kernel = m ake_kernel(model,q_vectors)189 kernel = model.make_kernel(q_vectors) 190 190 actual = call_kernel(kernel, pars) 191 191 else: 192 192 q_vectors = [np.array(x)] 193 kernel = m ake_kernel(model,q_vectors)193 kernel = model.make_kernel(q_vectors) 194 194 actual = call_kernel(kernel, pars) 195 195 -
sasmodels/models/cylinder.c
r26141cb r03cac08 3 3 double Iqxy(double qx, double qy, double sld, double solvent_sld, 4 4 double radius, double length, double theta, double phi); 5 6 #define INVALID(v) (v.radius<0 || v.length<0) 5 7 6 8 double form_volume(double radius, double length) … … 15 17 double length) 16 18 { 17 // TODO: return NaN if radius<0 or length<0?18 19 // precompute qr and qh to save time in the loop 19 20 const double qr = q*radius; … … 47 48 double phi) 48 49 { 49 // TODO: return NaN if radius<0 or length<0?50 50 double sn, cn; // slots to hold sincos function output 51 51 -
sasmodels/models/gel_fit.c
r30b4ddf r03cac08 1 double form_volume(void); 2 3 double Iq(double q, 4 double guinier_scale, 5 double lorentzian_scale, 6 double gyration_radius, 7 double fractal_exp, 8 double cor_length); 9 10 double Iqxy(double qx, double qy, 11 double guinier_scale, 12 double lorentzian_scale, 13 double gyration_radius, 14 double fractal_exp, 15 double cor_length); 16 17 static double _gel_fit_kernel(double q, 1 static double Iq(double q, 18 2 double guinier_scale, 19 3 double lorentzian_scale, … … 24 8 // Lorentzian Term 25 9 ////////////////////////double a(x[i]*x[i]*zeta*zeta); 26 double lorentzian_term = q*q*cor_length*cor_length;10 double lorentzian_term = square(q*cor_length); 27 11 lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 28 12 lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); … … 30 14 // Exponential Term 31 15 ////////////////////////double d(x[i]*x[i]*rg*rg); 32 double exp_term = q*q*gyration_radius*gyration_radius;16 double exp_term = square(q*gyration_radius); 33 17 exp_term = exp(-1.0*(exp_term/3.0)); 34 18 … … 37 21 return result; 38 22 } 39 double form_volume(void){40 // Unused, so free to return garbage.41 return NAN;42 }43 44 double Iq(double q,45 double guinier_scale,46 double lorentzian_scale,47 double gyration_radius,48 double fractal_exp,49 double cor_length)50 {51 return _gel_fit_kernel(q,52 guinier_scale,53 lorentzian_scale,54 gyration_radius,55 fractal_exp,56 cor_length);57 }58 59 // Iqxy is never called since no orientation or magnetic parameters.60 double Iqxy(double qx, double qy,61 double guinier_scale,62 double lorentzian_scale,63 double gyration_radius,64 double fractal_exp,65 double cor_length)66 {67 double q = sqrt(qx*qx + qy*qy);68 return _gel_fit_kernel(q,69 guinier_scale,70 lorentzian_scale,71 gyration_radius,72 fractal_exp,73 cor_length);74 }75 -
sasmodels/models/lib/Si.c
re7678b2 rba32cdd 1 // integral of sin(x)/x Taylor series approximated to w/i 0.1% 1 2 double Si(double x); 2 3 3 // integral of sin(x)/x Taylor series approximated to w/i 0.1%4 4 double Si(double x) 5 5 { -
sasmodels/models/lib/polevl.c
r0b05c24 rba32cdd 51 51 */ 52 52 53 double polevl( double x, constant double *coef, int N ) { 53 double polevl( double x, constant double *coef, int N ); 54 double p1evl( double x, constant double *coef, int N ); 55 56 57 double polevl( double x, constant double *coef, int N ) 58 { 54 59 55 60 int i = 0; … … 70 75 */ 71 76 72 double p1evl( double x, constant double *coef, int N ) { 77 double p1evl( double x, constant double *coef, int N ) 78 { 73 79 int i=0; 74 80 double ans = x+coef[i]; -
sasmodels/models/lib/sas_J0.c
ra776125 rba32cdd 49 49 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 50 50 */ 51 52 double sas_J0(double x); 51 53 52 54 /* Note: all coefficients satisfy the relative error criterion … … 177 179 }; 178 180 179 double sas_J0(double x) { 181 double sas_J0(double x) 182 { 180 183 181 184 //Cephes single precission -
sasmodels/models/lib/sas_J1.c
r19b9005 rba32cdd 39 39 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 40 40 */ 41 double sas_J1(double x); 42 double sas_J1c(double x); 41 43 42 44 constant double RPJ1[8] = { … … 135 137 }; 136 138 137 double sas_J1(double x) { 139 double sas_J1(double x) 140 { 138 141 139 142 //Cephes double pression function -
sasmodels/models/lib/sas_JN.c
ra776125 rba32cdd 47 47 Copyright 1984, 1987, 2000 by Stephen L. Moshier 48 48 */ 49 50 49 double sas_JN( int n, double x ); 51 50 52 double sas_JN( int n, double x ) { 51 double sas_JN( int n, double x ) 52 { 53 53 54 54 const double MACHEP = 1.11022302462515654042E-16; -
sasmodels/models/lib/sph_j1c.c
re6f1410 rba32cdd 7 7 * using double precision that are the source. 8 8 */ 9 double sph_j1c(double q); 9 10 10 11 // The choice of the number of terms in the series and the cutoff value for … … 43 44 #endif 44 45 45 double sph_j1c(double q);46 46 double sph_j1c(double q) 47 47 { -
sasmodels/models/lib/sphere_form.c
rad90df9 rba32cdd 1 inline double 2 sphere_volume(double radius) 1 double sphere_volume(double radius); 2 double sphere_form(double q, double radius, double sld, double solvent_sld); 3 4 double sphere_volume(double radius) 3 5 { 4 6 return M_4PI_3*cube(radius); 5 7 } 6 8 7 inline double 8 sphere_form(double q, double radius, double sld, double solvent_sld) 9 double sphere_form(double q, double radius, double sld, double solvent_sld) 9 10 { 10 11 const double fq = sphere_volume(radius) * sph_j1c(q*radius); -
sasmodels/models/lib/wrc_cyl.c
re7678b2 rba32cdd 2 2 Functions for WRC implementation of flexible cylinders 3 3 */ 4 double Sk_WR(double q, double L, double b); 5 6 4 7 static double 5 8 AlphaSquare(double x) … … 363 366 } 364 367 365 double Sk_WR(double q, double L, double b);366 368 double Sk_WR(double q, double L, double b) 367 369 { -
sasmodels/models/onion.c
rfdb1487 rce896fd 4 4 double thickness, double A) 5 5 { 6 const double vol = 4.0/3.0 * M_PI * r * r * r;6 const double vol = M_4PI_3 * cube(r); 7 7 const double qr = q * r; 8 8 const double alpha = A * r/thickness; … … 19 19 double sinqr, cosqr; 20 20 SINCOS(qr, sinqr, cosqr); 21 fun = -3.0 (21 fun = -3.0*( 22 22 ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 23 23 - (alpha*sinqr/qr - cosqr) … … 32 32 f_linear(double q, double r, double sld, double slope) 33 33 { 34 const double vol = 4.0/3.0 * M_PI * r * r * r;34 const double vol = M_4PI_3 * cube(r); 35 35 const double qr = q * r; 36 36 const double bes = sph_j1c(qr); … … 52 52 { 53 53 const double bes = sph_j1c(q * r); 54 const double vol = 4.0/3.0 * M_PI * r * r * r;54 const double vol = M_4PI_3 * cube(r); 55 55 return sld * vol * bes; 56 56 } … … 64 64 r += thickness[i]; 65 65 } 66 return 4.0/3.0 * M_PI * r * r * r;66 return M_4PI_3*cube(r); 67 67 } 68 68 69 69 static double 70 Iq(double q, double core_sld, double core_radius, double solvent_sld,71 double n, double in_sld[], double out_sld[], double thickness[],70 Iq(double q, double sld_core, double core_radius, double sld_solvent, 71 double n, double sld_in[], double sld_out[], double thickness[], 72 72 double A[]) 73 73 { 74 74 int i; 75 r = core_radius;76 f = f_constant(q, r, core_sld);75 double r = core_radius; 76 double f = f_constant(q, r, sld_core); 77 77 for (i=0; i<n; i++){ 78 78 const double r0 = r; … … 92 92 } 93 93 } 94 f -= f_constant(q, r, s olvent_sld);95 f2 = f * f * 1.0e-4;94 f -= f_constant(q, r, sld_solvent); 95 const double f2 = f * f * 1.0e-4; 96 96 97 97 return f2; -
sasmodels/models/onion.py
rb0696e1 rb0696e1 293 293 294 294 # ["name", "units", default, [lower, upper], "type","description"], 295 parameters = [[" core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",295 parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 296 296 "Core scattering length density"], 297 297 ["core_radius", "Ang", 200., [0, inf], "volume", 298 298 "Radius of the core"], 299 ["s olvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "",299 ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 300 300 "Solvent scattering length density"], 301 301 ["n", "", 1, [0, 10], "volume", 302 302 "number of shells"], 303 [" in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",303 ["sld_in[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 304 304 "scattering length density at the inner radius of shell k"], 305 [" out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",305 ["sld_out[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 306 306 "scattering length density at the outer radius of shell k"], 307 307 ["thickness[n]", "Ang", 40., [0, inf], "volume", … … 311 311 ] 312 312 313 #source = ["lib/sph_j1c.c", "onion.c"] 314 315 def Iq(q, *args, **kw): 316 return q 317 318 def Iqxy(qx, *args, **kw): 319 return qx 320 321 322 def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 313 source = ["lib/sph_j1c.c", "onion.c"] 314 315 #def Iq(q, *args, **kw): 316 # return q 317 318 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 319 def profile(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 323 320 """ 324 321 SLD profile … … 374 371 375 372 demo = { 376 "s olvent_sld": 2.2,377 " core_sld": 1.0,373 "sld_solvent": 2.2, 374 "sld_core": 1.0, 378 375 "core_radius": 100, 379 376 "n": 4, 380 " in_sld": [0.5, 1.5, 0.9, 2.0],381 " out_sld": [nan, 0.9, 1.2, 1.6],377 "sld_in": [0.5, 1.5, 0.9, 2.0], 378 "sld_out": [nan, 0.9, 1.2, 1.6], 382 379 "thickness": [50, 75, 150, 75], 383 380 "A": [0, -1, 1e-4, 1], -
sasmodels/models/rpa.py
raad336c r69aa451 86 86 # ["name", "units", default, [lower, upper], "type","description"], 87 87 parameters = [ 88 ["case_num", CASES, 0, [0, 10], "", "Component organization"],88 ["case_num", "", 1, CASES, "", "Component organization"], 89 89 90 ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], 91 ["Phia", "", 0.25, [0, 1], "", "volume fraction"], 92 ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"], 93 ["La", "fm", 10.0, [-inf, inf], "", "scattering length"], 94 ["ba", "Ang", 5.0, [0, inf], "", "segment length"], 95 96 ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"], 97 ["Phib", "", 0.25, [0, 1], "", "volume fraction"], 98 ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"], 99 ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"], 100 ["bb", "Ang", 5.0, [0, inf], "", "segment length"], 101 102 ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"], 103 ["Phic", "", 0.25, [0, 1], "", "volume fraction"], 104 ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"], 105 ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"], 106 ["bc", "Ang", 5.0, [0, inf], "", "segment length"], 107 108 ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"], 109 ["Phid", "", 0.25, [0, 1], "", "volume fraction"], 110 ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"], 111 ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"], 112 ["bd", "Ang", 5.0, [0, inf], "", "segment length"], 90 ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 91 ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 92 ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 93 ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 94 ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 113 95 114 96 ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"], -
sasmodels/models/spherical_sld.py
r6f0e04f r6f0e04f 163 163 title = "Sperical SLD intensity calculation" 164 164 description = """ 165 166 167 165 I(q) = 166 background = Incoherent background [1/cm] 167 """ 168 168 category = "sphere-based" 169 170 INTERFACE_CASES = ["Erf", "RPower", "LPower", "RExp", "LExp"] 169 171 170 172 # pylint: disable=bad-whitespace, line-too-long 171 173 # ["name", "units", default, [lower, upper], "type", "description"], 172 parameters = [["n_shells", "", 1, [0, 9],"", "number of shells"],173 ["radius_core", "Ang", 50.0, [0, inf],"", "intern layer thickness"],174 ["sld_core", "1e-6/Ang^2", 2.07, [-inf, inf],"", "sld function flat"],175 ["sld_ flat[n]", "1e-6/Ang^2", 4.06, [-inf, inf], "", "sld function flat"],176 [" thick_flat[n]", "Ang", 100.0, [0, inf], "", "flat layer_thickness"],177 [" func_inter[n]", "", 0, [0, 4], "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"],178 [" thick_inter[n]", "Ang", 50.0, [0, inf], "", "intern layer thickness"],179 [" inter_nu[n]", "", 2.5, [-inf, inf],"", "steepness parameter"],180 [" npts_inter", "", 35, [0, 35], "", "number of points in each sublayer Must be odd number"],181 [" sld_solvent", "1e-6/Ang^2", 1.0, [-inf, inf], "", "sld function solvent"],174 parameters = [["n_shells", "", 1, [0, 9], "", "number of shells"], 175 ["radius_core", "Ang", 50.0, [0, inf], "", "intern layer thickness"], 176 ["sld_core", "1e-6/Ang^2", 2.07, [-inf, inf], "", "sld function flat"], 177 ["sld_solvent", "1e-6/Ang^2", 1.0, [-inf, inf], "", "sld function solvent"], 178 ["sld_flat[n_shells]", "1e-6/Ang^2", 4.06, [-inf, inf], "", "sld function flat"], 179 ["thick_flat[n_shells]", "Ang", 100.0, [0, inf], "", "flat layer_thickness"], 180 ["func_inter[n_shells]", "", 0, INTERFACE_CASES, "", "shape of the interface between shells"], 181 ["nu_inter[n_shells]", "", 2.5, [-inf, inf], "", "steepness parameter"], 182 ["thick_inter[n_shells]", "Ang", 50.0, [0, inf], "", "intern layer thickness"], 183 ["n_points_inter", "", 35, [0, 35], "", "number of points in each sublayer Must be odd number"], 182 184 ] 183 185 # pylint: enable=bad-whitespace, line-too-long 186 184 187 #source = ["lib/librefl.c", "lib/sph_j1c.c", "spherical_sld.c"] 185 186 188 def Iq(q, *args, **kw): 187 189 return q 188 189 def Iqxy(qx, *args, **kw):190 return qx191 192 190 193 191 demo = dict( … … 196 194 solvent_sld=1.0, 197 195 background=0.0, 198 n pts_inter=35.0,196 n_points_inter=35.0, 199 197 ) 200 198 … … 210 208 #func_inter='func_inter', 211 209 #thick_inter='thick_inter', 212 # inter_nu='nu_inter',213 # npts_inter='npts_inter',210 #nu_inter='nu_inter', 211 #points_inter='npts_inter', 214 212 sld_solvent='sld_solv', 215 213 ) … … 218 216 tests = [ 219 217 # Accuracy tests based on content in test/utest_extra_models.py 220 [{'n pts_iter':35,221 222 223 224 225 226 227 228 229 230 231 232 233 }, 0.001, 0.001],218 [{'n_points_inter':35, 219 'sld_solv':1, 220 'radius_core':50.0, 221 'sld_core':2.07, 222 'func_inter2':0.0, 223 'thick_inter2':50, 224 'nu_inter2':2.5, 225 'sld_flat2':4, 226 'thick_flat2':100, 227 'func_inter1':0.0, 228 'thick_inter1':50, 229 'nu_inter1':2.5, 230 'background': 0.0, 231 }, 0.001, 0.001], 234 232 ] -
sasmodels/product.py
r3c6d5bc r3c6d5bc 14 14 15 15 from .core import call_ER_VR 16 from .generate import process_parameters17 16 18 17 SCALE=0 … … 66 65 model_info['oldpars'] = oldpars 67 66 model_info['composition'] = ('product', [p_info, s_info]) 68 process_parameters(model_info)69 67 return model_info 70 68 … … 101 99 # a parameter map. 102 100 par_map = {} 103 p_info = p_kernel.info['par type']104 s_info = s_kernel.info['par type']101 p_info = p_kernel.info['par_type'] 102 s_info = s_kernel.info['par_type'] 105 103 vol_pars = set(p_info['volume']) 106 104 if dim == '2d': -
sasmodels/resolution.py
r486fcf6 r486fcf6 477 477 """ 478 478 from sasmodels import core 479 kernel = core.make_kernel(form,[q])479 kernel = form.make_kernel([q]) 480 480 theory = core.call_kernel(kernel, pars) 481 481 kernel.release() … … 502 502 from scipy.integrate import romberg 503 503 504 if any(k not in form.info['defaults'] for k in pars.keys()):505 keys = set(form.info['defaults'].keys())506 extra = set(pars.keys()) - keys504 par_set = set([p.name for p in form.info['parameters']]) 505 if any(k not in par_set for k in pars.keys()): 506 extra = set(pars.keys()) - par_set 507 507 raise ValueError("bad parameters: [%s] not in [%s]"% 508 508 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) … … 556 556 from scipy.integrate import romberg 557 557 558 if any(k not in form.info['defaults'] for k in pars.keys()):559 keys = set(form.info['defaults'].keys())560 extra = set(pars.keys()) - keys558 par_set = set([p.name for p in form.info['parameters']]) 559 if any(k not in par_set for k in pars.keys()): 560 extra = set(pars.keys()) - par_set 561 561 raise ValueError("bad parameters: [%s] not in [%s]"% 562 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 562 (", ".join(sorted(extra)), 563 ", ".join(sorted(pars.keys())))) 563 564 564 565 _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) … … 693 694 def _eval_sphere(self, pars, resolution): 694 695 from sasmodels import core 695 kernel = core.make_kernel(self.model,[resolution.q_calc])696 kernel = self.model.make_kernel([resolution.q_calc]) 696 697 theory = core.call_kernel(kernel, pars) 697 698 result = resolution.apply(theory) … … 1062 1063 model = core.build_model(model_info) 1063 1064 1064 kernel = core.make_kernel(model,[resolution.q_calc])1065 kernel = model.make_kernel([resolution.q_calc]) 1065 1066 theory = core.call_kernel(kernel, pars) 1066 1067 Iq = resolution.apply(theory) -
sasmodels/sasview_model.py
r787be86 rfb5914f 21 21 22 22 from . import core 23 from . import generate23 from . import weights 24 24 25 25 def standard_models(): … … 32 32 33 33 Returns a class that can be used directly as a sasview model.t 34 35 Defaults to using the new name for a model. Setting36 *namestyle='oldname'* will produce a class with a name37 compatible with SasView.38 34 """ 39 35 model_info = core.load_model_info(model_name) … … 56 52 """ 57 53 def __init__(self): 58 self._ kernel = None54 self._model = None 59 55 model_info = self._model_info 56 parameters = model_info['parameters'] 60 57 61 58 self.name = model_info['name'] 62 self.oldname = model_info['oldname']63 59 self.description = model_info['description'] 64 60 self.category = None 65 self.multiplicity_info = None 66 self.is_multifunc = False 61 #self.is_multifunc = False 62 for p in parameters.kernel_parameters: 63 if p.is_control: 64 profile_axes = model_info['profile_axes'] 65 self.multiplicity_info = [ 66 p.limits[1], p.name, p.choices, profile_axes[0] 67 ] 68 break 69 else: 70 self.multiplicity_info = [] 67 71 68 72 ## interpret the parameters … … 71 75 self.params = collections.OrderedDict() 72 76 self.dispersion = dict() 73 partype = model_info['partype'] 74 75 for p in model_info['parameters']: 77 78 self.orientation_params = [] 79 self.magnetic_params = [] 80 self.fixed = [] 81 for p in parameters.user_parameters(): 76 82 self.params[p.name] = p.default 77 83 self.details[p.name] = [p.units] + p.limits 78 79 for name in partype['pd-2d']: 80 self.dispersion[name] = { 81 'width': 0, 82 'npts': 35, 83 'nsigmas': 3, 84 'type': 'gaussian', 85 } 86 87 self.orientation_params = ( 88 partype['orientation'] 89 + [n + '.width' for n in partype['orientation']] 90 + partype['magnetic']) 91 self.magnetic_params = partype['magnetic'] 92 self.fixed = [n + '.width' for n in partype['pd-2d']] 84 if p.polydisperse: 85 self.dispersion[p.name] = { 86 'width': 0, 87 'npts': 35, 88 'nsigmas': 3, 89 'type': 'gaussian', 90 } 91 if p.type == 'orientation': 92 self.orientation_params.append(p.name) 93 self.orientation_params.append(p.name+".width") 94 self.fixed.append(p.name+".width") 95 if p.type == 'magnetic': 96 self.orientation_params.append(p.name) 97 self.magnetic_params.append(p.name) 98 self.fixed.append(p.name+".width") 99 93 100 self.non_fittable = [] 94 101 … … 109 116 def __get_state__(self): 110 117 state = self.__dict__.copy() 111 model_id = self._model_info['id']112 118 state.pop('_kernel') 113 119 # May need to reload model info on set state since it has pointers … … 118 124 def __set_state__(self, state): 119 125 self.__dict__ = state 120 self._ kernel = None126 self._model = None 121 127 122 128 def __str__(self): … … 207 213 def getDispParamList(self): 208 214 """ 209 Return a list of all availableparameters for the model215 Return a list of polydispersity parameters for the model 210 216 """ 211 217 # TODO: fix test so that parameter order doesn't matter 212 ret = ['%s.%s' % (d.lower(), p) 213 for d in self._model_info['partype']['pd-2d'] 214 for p in ('npts', 'nsigmas', 'width')] 218 ret = ['%s.%s' % (p.name.lower(), ext) 219 for p in self._model_info['parameters'].user_parameters() 220 for ext in ('npts', 'nsigmas', 'width') 221 if p.polydisperse] 215 222 #print(ret) 216 223 return ret … … 285 292 # Check whether we have a list of ndarrays [qx,qy] 286 293 qx, qy = qdist 287 partype = self._model_info['partype'] 288 if not partype['orientation'] and not partype['magnetic']: 294 if not self._model_info['parameters'].has_2d: 289 295 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 290 296 else: … … 308 314 to the card for each evaluation. 309 315 """ 310 if self._ kernel is None:311 self._ kernel = core.build_model(self._model_info)316 if self._model is None: 317 self._model = core.build_model(self._model_info, platform='dll') 312 318 q_vectors = [np.asarray(q) for q in args] 313 fn = self._kernel(q_vectors) 314 pars = [self.params[v] for v in fn.fixed_pars] 315 pd_pars = [self._get_weights(p) for p in fn.pd_pars] 316 result = fn(pars, pd_pars, self.cutoff) 317 fn.q_input.release() 318 fn.release() 319 kernel = self._model.make_kernel(q_vectors) 320 pairs = [self._get_weights(p) 321 for p in self._model_info['parameters'].call_parameters] 322 details, weights, values = core.build_details(kernel, pairs) 323 return kernel(details, weights, values, cutoff=self.cutoff) 324 kernel.q_input.release() 325 kernel.release() 319 326 return result 320 327 … … 389 396 def _get_weights(self, par): 390 397 """ 391 Return dispersion weights 392 :param par parameter name 393 """ 394 from . import weights 395 396 relative = self._model_info['partype']['pd-rel'] 397 limits = self._model_info['limits'] 398 dis = self.dispersion[par] 399 value, weight = weights.get_weights( 400 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 401 self.params[par], limits[par], par in relative) 402 return value, weight / np.sum(weight) 403 398 Return dispersion weights for parameter 399 """ 400 if par.polydisperse: 401 dis = self.dispersion[par.name] 402 value, weight = weights.get_weights( 403 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 404 self.params[par.name], par.limits, par.relative_pd) 405 return value, weight / np.sum(weight) 406 else: 407 return [self.params[par.name]], [] 408 409 def test_model(): 410 Cylinder = make_class('cylinder') 411 cylinder = Cylinder() 412 return cylinder.evalDistribution([0.1,0.1]) 413 414 if __name__ == "__main__": 415 print("cylinder(0.1,0.1)=%g"%test_model())
Note: See TracChangeset
for help on using the changeset viewer.