Changes in / [34edbb8:e9b1663d] in sasmodels
- Files:
-
- 6 added
- 16 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
r9217ef8 r9217ef8 81 81 from bumps.names import Parameter 82 82 83 pars = {} 83 pars = {} # => floating point parameters 84 84 for p in model_info['parameters']: 85 85 value = kwargs.pop(p.name, p.default) 86 86 pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 87 for name in model_info['par type']['pd-2d']:87 for name in model_info['par_type']['pd']: 88 88 for xpart, xdefault, xlimits in [ 89 89 ('_pd', 0., pars[name].limits), … … 95 95 pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 96 96 97 pd_types = {} 98 for name in model_info['par type']['pd-2d']:97 pd_types = {} # => distribution names 98 for name in model_info['par_type']['pd']: 99 99 xname = name + '_pd_type' 100 100 xvalue = kwargs.pop(xname, 'gaussian') -
sasmodels/compare.py
r1671636 r69aa451 309 309 exclude = lambda n: False 310 310 else: 311 partype = model_info['partype'] 312 par1d = set(partype['fixed-1d']+partype['pd-1d']) 311 par1d = model_info['par_type']['1d'] 313 312 exclude = lambda n: n not in par1d 314 313 lines = [] … … 454 453 """ 455 454 # initialize the code so time is more accurate 456 value = calculator(**suppress_pd(pars)) 455 if Nevals > 1: 456 value = calculator(**suppress_pd(pars)) 457 457 toc = tic() 458 458 for _ in range(max(Nevals, 1)): # make sure there is at least one eval … … 877 877 pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 878 878 if not opts['is2d']: 879 active = [base + ext 880 for base in model_info['partype']['pd-1d'] 881 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 882 active.extend(model_info['partype']['fixed-1d']) 883 for k in active: 884 v = pars[k] 885 v.range(*parameter_range(k, v.value)) 879 for p in model_info['parameters'].type['1d']: 880 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 881 k = p.name+ext 882 v = pars.get(k, None) 883 if v is not None: 884 v.range(*parameter_range(k, v.value)) 886 885 else: 887 886 for k, v in pars.items(): -
sasmodels/core.py
re79f0a5 r69aa451 26 26 __all__ = [ 27 27 "list_models", "load_model_info", "precompile_dll", 28 "build_model", " make_kernel", "call_kernel", "call_ER_VR",28 "build_model", "call_kernel", "call_ER_VR", 29 29 ] 30 30 … … 167 167 168 168 169 def make_kernel(model, q_vectors): 170 """ 171 Return a computation kernel from the model definition and the q input. 172 """ 173 return model(q_vectors) 174 175 def get_weights(model_info, pars, name): 169 def get_weights(parameter, values): 176 170 """ 177 171 Generate the distribution for parameter *name* given the parameter values … … 181 175 from the *pars* dictionary for parameter value and parameter dispersion. 182 176 """ 183 relative = name in model_info['partype']['pd-rel'] 184 limits = model_info['limits'][name] 185 disperser = pars.get(name+'_pd_type', 'gaussian') 186 value = pars.get(name, model_info['defaults'][name]) 187 npts = pars.get(name+'_pd_n', 0) 188 width = pars.get(name+'_pd', 0.0) 189 nsigma = pars.get(name+'_pd_nsigma', 3.0) 177 value = values.get(parameter.name, parameter.default) 178 if parameter.type not in ('volume', 'orientation'): 179 return np.array([value]), np.array([1.0]) 180 relative = parameter.type == 'volume' 181 limits = parameter.limits 182 disperser = values.get(parameter.name+'_pd_type', 'gaussian') 183 npts = values.get(parameter.name+'_pd_n', 0) 184 width = values.get(parameter.name+'_pd', 0.0) 185 nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 190 186 value, weight = weights.get_weights( 191 187 disperser, npts, width, nsigma, value, limits, relative) … … 208 204 def call_kernel(kernel, pars, cutoff=0, mono=False): 209 205 """ 210 Call *kernel* returned from :func:`make_kernel`with parameters *pars*.206 Call *kernel* returned from *model.make_kernel* with parameters *pars*. 211 207 212 208 *cutoff* is the limiting value for the product of dispersion weights used … … 216 212 with an error of about 1%, which is usually less than the measurement 217 213 uncertainty. 218 """ 219 fixed_pars = [pars.get(name, kernel.info['defaults'][name])220 for name in kernel.fixed_pars]214 215 *mono* is True if polydispersity should be set to none on all parameters. 216 """ 221 217 if mono: 222 pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 223 for name in kernel.pd_pars] 224 else: 225 pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 226 return kernel(fixed_pars, pd_pars, cutoff=cutoff) 218 active = lambda name: False 219 elif kernel.dim == '1d': 220 pars_1d = set(p.name for p in kernel.info['parameters'].type['1d']) 221 active = lambda name: name in pars_1d 222 elif kernel.dim == '2d': 223 pars_2d = set(p.name for p in kernel.info['parameters'].type['2d']) 224 active = lambda name: name in pars_1d 225 else: 226 active = lambda name: True 227 228 vw_pairs = [(get_weights(p, pars) if active(p.name) else ([p.default], [1])) 229 for p in kernel.info['parameters']] 230 values, weights = zip(*vw_pairs) 231 232 if max([len(w) for w in weights]) > 1: 233 details = generate.poly_details(kernel.info, weights) 234 else: 235 details = kernel.info['mono_details'] 236 237 weights, values = [np.hstack(v) for v in (weights, values)] 238 239 weights = weights.astype(dtype=kernel.dtype) 240 values = values.astype(dtype=kernel.dtype) 241 return kernel(details, weights, values, cutoff) 227 242 228 243 def call_ER_VR(model_info, vol_pars): … … 247 262 248 263 249 def call_ER( info, pars):250 """ 251 Call the model ER function using * pars*.252 * info* is either *model.info* if you have a loaded model, or *kernel.info*253 if youhave a model kernel prepared for evaluation.254 """ 255 ER = info.get('ER', None)264 def call_ER(model_info, values): 265 """ 266 Call the model ER function using *values*. *model_info* is either 267 *model.info* if you have a loaded model, or *kernel.info* if you 268 have a model kernel prepared for evaluation. 269 """ 270 ER = model_info.get('ER', None) 256 271 if ER is None: 257 272 return 1.0 258 273 else: 259 vol_pars = [get_weights(info, pars, name) 260 for name in info['partype']['volume']] 274 vol_pars = [get_weights(parameter, values) 275 for parameter in model_info['parameters'] 276 if parameter.type == 'volume'] 261 277 value, weight = dispersion_mesh(vol_pars) 262 278 individual_radii = ER(*value) … … 264 280 return np.sum(weight*individual_radii) / np.sum(weight) 265 281 266 def call_VR( info, pars):282 def call_VR(model_info, values): 267 283 """ 268 284 Call the model VR function using *pars*. … … 270 286 if you have a model kernel prepared for evaluation. 271 287 """ 272 VR = info.get('VR', None)288 VR = model_info.get('VR', None) 273 289 if VR is None: 274 290 return 1.0 275 291 else: 276 vol_pars = [get_weights(info, pars, name) 277 for name in info['partype']['volume']] 292 vol_pars = [get_weights(parameter, values) 293 for parameter in model_info['parameters'] 294 if parameter.type == 'volume'] 278 295 value, weight = dispersion_mesh(vol_pars) 279 296 whole, part = VR(*value) -
sasmodels/direct_model.py
r02e70ff r48fbd50 25 25 import numpy as np 26 26 27 from .core import make_kernel28 27 from .core import call_kernel, call_ER_VR 29 28 from . import sesans … … 66 65 self.data_type = 'Iq' 67 66 68 partype = model.info['partype']69 70 67 if self.data_type == 'sesans': 71 68 … … 81 78 q_mono = sesans.make_all_q(data) 82 79 elif self.data_type == 'Iqxy': 80 partype = model.info['par_type'] 83 81 if not partype['orientation'] and not partype['magnetic']: 84 82 raise ValueError("not 2D without orientation or magnetic parameters") … … 154 152 def _calc_theory(self, pars, cutoff=0.0): 155 153 if self._kernel is None: 156 self._kernel = make_kernel(self._model, self._kernel_inputs) # pylint: disable=attribute-dedata_type 157 self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None 154 self._kernel = self._model.make_kernel(self._kernel_inputs) 155 self._kernel_mono = ( 156 self._model.make_kernel(self._kernel_mono_inputs) 157 if self._kernel_mono_inputs else None 158 ) 158 159 159 160 Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 160 Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None 161 Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 162 if self._kernel_mono_inputs else None) 161 163 if self.data_type == 'sesans': 162 164 result = sesans.transform(self._data, -
sasmodels/generate.py
r9ef9dd9 r69aa451 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. 183 #TODO: determine which functions are useful outside of generate 184 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 215 185 216 186 import sys 217 187 from os.path import abspath, dirname, join as joinpath, exists, basename, \ 218 splitext 188 splitext, getmtime 219 189 import re 220 190 import string 221 191 import warnings 222 from collections import namedtuple223 192 224 193 import numpy as np 225 194 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') 195 from .modelinfo import ModelInfo, Parameter, make_parameter_table 196 197 # TODO: identify model files which have changed since loading and reload them. 198 199 TEMPLATE_ROOT = dirname(__file__) 200 201 MAX_PD = 4 233 202 234 203 F16 = np.dtype('float16') … … 239 208 except TypeError: 240 209 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 210 248 211 # Conversion from units defined in the parameter table for each model … … 338 301 raise ValueError("%r not found in %s" % (filename, search_path)) 339 302 303 340 304 def model_sources(model_info): 341 305 """ … … 346 310 return [_search(search_path, f) for f in model_info['source']] 347 311 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 """ 312 def timestamp(model_info): 313 """ 314 Return a timestamp for the model corresponding to the most recently 315 changed file or dependency. 316 """ 317 source_files = (model_sources(model_info) 318 + model_templates() 319 + [model_info['filename']]) 320 newest = max(getmtime(f) for f in source_files) 321 return newest 361 322 362 323 def convert_type(source, dtype): … … 369 330 if dtype == F16: 370 331 fbytes = 2 371 source = _ F16_PRAGMA + _convert_type(source, "half", "f")332 source = _convert_type(source, "float", "f") 372 333 elif dtype == F32: 373 334 fbytes = 4 … … 375 336 elif dtype == F64: 376 337 fbytes = 8 377 source = _F64_PRAGMA + source # Sourceis already double338 # no need to convert if it is already double 378 339 elif dtype == F128: 379 340 fbytes = 16 … … 418 379 419 380 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];\ 381 _template_cache = {} 382 def load_template(filename): 383 path = joinpath(TEMPLATE_ROOT, filename) 384 mtime = getmtime(path) 385 if filename not in _template_cache or mtime > _template_cache[filename][0]: 386 with open(path) as fid: 387 _template_cache[filename] = (mtime, fid.read(), path) 388 return _template_cache[filename][1] 389 390 def model_templates(): 391 # TODO: fails DRY; templates are listed in two places. 392 # should instead have model_info contain a list of paths 393 return [joinpath(TEMPLATE_ROOT, filename) 394 for filename in ('kernel_header.c', 'kernel_iq.c')] 395 396 397 _FN_TEMPLATE = """\ 398 double %(name)s(%(pars)s); 399 double %(name)s(%(pars)s) { 400 %(body)s 401 } 402 403 424 404 """ 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 405 406 def _gen_fn(name, pars, body): 407 """ 408 Generate a function given pars and body. 409 410 Returns the following string:: 411 412 double fn(double a, double b, ...); 413 double fn(double a, double b, ...) { 414 .... 415 } 416 """ 417 par_decl = ', '.join(p.as_argument() for p in pars) if pars else 'void' 418 return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 419 420 def _call_pars(prefix, pars): 421 """ 422 Return a list of *prefix.parameter* from parameter items. 423 """ 424 return [p.as_call_reference(prefix) for p in pars] 425 426 _IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 427 flags=re.MULTILINE) 428 def _have_Iqxy(sources): 429 """ 430 Return true if any file defines Iqxy. 431 432 Note this is not a C parser, and so can be easily confused by 433 non-standard syntax. Also, it will incorrectly identify the following 434 as having Iqxy:: 435 436 /* 437 double Iqxy(qx, qy, ...) { ... fill this in later ... } 438 */ 439 440 If you want to comment out an Iqxy function, use // on the front of the 441 line instead. 442 """ 443 for code in sources: 444 if _IQXY_PATTERN.search(code): 445 return True 446 else: 447 return False 448 444 449 def make_source(model_info): 445 450 """ … … 460 465 # for computing volume even if we allow non-disperse volume parameters. 461 466 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 467 partable = model_info['parameters'] 468 469 # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 470 # Note that scale and volume are not possible types. 471 472 # Load templates and user code 473 kernel_header = load_template('kernel_header.c') 474 kernel_code = load_template('kernel_iq.c') 475 user_code = [open(f).read() for f in model_sources(model_info)] 476 477 # Build initial sources 478 source = [kernel_header] + user_code 479 480 vol_parameters = partable.kernel_pars('volume') 481 iq_parameters = partable.kernel_pars('1d') 482 iqxy_parameters = partable.kernel_pars('2d') 483 484 # Make parameters for q, qx, qy so that we can use them in declarations 485 q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 486 # Generate form_volume function, etc. from body only 497 487 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'))) 488 pars = vol_parameters 489 source.append(_gen_fn('form_volume', pars, model_info['form_volume'])) 530 490 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'))) 491 pars = [q] + iq_parameters 492 source.append(_gen_fn('Iq', pars, model_info['Iq'])) 559 493 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 } 494 pars = [qx, qy] + iqxy_parameters 495 source.append(_gen_fn('Iqxy', pars, model_info['Iqxy'])) 496 497 # Define the parameter table 498 source.append("#define PARAMETER_TABLE \\") 499 source.append("\\\n".join(p.as_definition() 500 for p in model_info['parameters'][2:])) 501 502 # Define the function calls 503 if vol_parameters: 504 refs = _call_pars("v.", vol_parameters) 505 call_volume = "#define CALL_VOLUME(v) form_volume(%s)" % (",".join(refs)) 506 else: 507 # Model doesn't have volume. We could make the kernel run a little 508 # faster by not using/transferring the volume normalizations, but 509 # the ifdef's reduce readability more than is worthwhile. 510 call_volume = "#define CALL_VOLUME(v) 0.0" 511 source.append(call_volume) 512 513 refs = ["q[i]"] + _call_pars("v.", iq_parameters) 514 call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs)) 515 if _have_Iqxy(user_code): 516 # Call 2D model 517 refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v.", iqxy_parameters) 518 call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs)) 519 else: 520 # Call 1D model with sqrt(qx^2 + qy^2) 521 warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 522 # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 523 pars_sqrt = ["sqrt(q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:] 524 call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt)) 525 526 # Fill in definitions for numbers of parameters 527 source.append("#define MAX_PD %s"%model_info['max_pd']) 528 source.append("#define NPARS %d"%(len(partable.kernel_pars()))) 529 530 # TODO: allow mixed python/opencl kernels? 531 532 # define the Iq kernel 533 source.append("#define KERNEL_NAME %s_Iq"%model_info['name']) 534 source.append(call_iq) 535 source.append(kernel_code) 536 source.append("#undef CALL_IQ") 537 source.append("#undef KERNEL_NAME") 538 539 # define the Iqxy kernel from the same source with different #defines 540 source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name']) 541 source.append(call_iqxy) 542 source.append(kernel_code) 543 source.append("#undef CALL_IQ") 544 source.append("#undef KERNEL_NAME") 545 546 return '\n'.join(source) 583 547 584 548 def categorize_parameters(pars): 585 549 """ 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 550 Categorize the parameters by use: 551 552 * *pd* list of polydisperse parameters in order; gui should test whether 553 they are in *2d* or *magnetic* as appropriate for the data 554 * *1d* set of parameters that are used to compute 1D patterns 555 * *2d* set of parameters that are used to compute 2D patterns (which 556 includes all 1D parameters) 557 * *magnetic* set of parameters that are used to compute magnetic 558 patterns (which includes all 1D and 2D parameters) 559 * *pd_relative* is the set of parameters with relative distribution 560 width (e.g., radius +/- 10%) rather than absolute distribution 561 width (e.g., theta +/- 6 degrees). 562 * *theta_par* is the index of the polar angle polydispersion parameter 563 or -1 if no such parameter exists 564 """ 565 par_set = {} 632 566 633 567 def process_parameters(model_info): … … 635 569 Process parameter block, precalculating parameter details. 636 570 """ 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) 571 partable = model_info['parameters'] 651 572 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'] 573 model_info['demo'] = partable.defaults 574 575 # Don't use more polydisperse parameters than are available in the model 576 # Note: we can do polydispersity on arbitrary parameters, so it is not 577 # clear that this is a good idea; it does however make the poly_details 578 # code easier to write, so we will leave it in for now. 579 model_info['max_pd'] = min(partable.num_pd, MAX_PD) 580 581 def mono_details(model_info): 582 # TODO: move max_pd into ParameterTable? 583 max_pd = model_info['max_pd'] 584 pars = model_info['parameters'].kernel_pars() 585 npars = len(pars) 586 par_offset = 5*max_pd 587 constants_offset = par_offset + 3*npars 588 589 details = np.zeros(constants_offset + 2, 'int32') 590 details[0*max_pd:1*max_pd] = range(max_pd) # pd_par: arbitrary order; use first 591 details[1*max_pd:2*max_pd] = [1]*max_pd # pd_length: only one element 592 details[2*max_pd:3*max_pd] = range(max_pd) # pd_offset: consecutive 1.0 weights 593 details[3*max_pd:4*max_pd] = [1]*max_pd # pd_stride: vectors of length 1 594 details[4*max_pd:5*max_pd] = [0]*max_pd # pd_isvol: doens't matter if no norm 595 details[par_offset+0*npars:par_offset+1*npars] = range(2, npars+2) # par_offset: skip scale and background 596 details[par_offset+1*npars:par_offset+2*npars] = [0]*npars # no coordination 597 #details[p+npars] = 1 # par_coord[0] is coordinated with the first par? 598 details[par_offset+2*npars:par_offset+3*npars] = 0 # fast coord with 0 599 details[constants_offset] = 1 # fast_coord_count: one fast index 600 details[constants_offset+1] = -1 # theta_par: None 601 return details 602 603 def poly_details(model_info, weights): 604 weights = weights[2:] 605 606 # TODO: move max_pd into ParameterTable? 607 max_pd = model_info['max_pd'] 608 pars = model_info['parameters'].kernel_pars 609 npars = len(pars) 610 par_offset = 5*max_pd 611 constants_offset = par_offset + 3*npars 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 print (pd_length) 617 print (weights) 618 pd_offset = np.cumsum(np.hstack((0, pd_length))) 619 pd_isvol = np.array([p.type=='volume' for p in pars]) 620 idx = np.argsort(pd_length)[::-1][:max_pd] 621 print (idx) 622 pd_stride = np.cumprod(np.hstack((1, pd_length[idx][:-1]))) 623 par_offsets = np.cumsum(np.hstack((2, pd_length)))[:-1] 624 625 theta_par = -1 626 if 'theta_par' in model_info: 627 theta_par = model_info['theta_par'] 628 if theta_par >= 0 and pd_length[theta_par] <= 1: 629 theta_par = -1 630 631 details = np.empty(constants_offset + 2, 'int32') 632 details[0*max_pd:1*max_pd] = idx # pd_par 633 details[1*max_pd:2*max_pd] = pd_length[idx] 634 details[2*max_pd:3*max_pd] = pd_offset[idx] 635 details[3*max_pd:4*max_pd] = pd_stride 636 details[4*max_pd:5*max_pd] = pd_isvol[idx] 637 details[par_offset+0*npars:par_offset+1*npars] = par_offsets 638 details[par_offset+1*npars:par_offset+2*npars] = 0 # no coordination for most 639 details[par_offset+2*npars:par_offset+3*npars] = 0 # no fast coord with 0 640 coord_offset = par_offset+1*npars 641 for k,parameter_num in enumerate(idx): 642 details[coord_offset+parameter_num] = 2**k 643 details[constants_offset] = 1 # fast_coord_count: one fast index 644 details[constants_offset+1] = theta_par 645 print ("details",details) 646 return details 647 648 def constrained_poly_details(model_info, weights, constraints): 649 # Need to find the independently varying pars and sort them 650 # Need to build a coordination list for the dependent variables 651 # Need to generate a constraints function which takes values 652 # and weights, returning par blocks 653 raise NotImplementedError("Can't handle constraints yet") 654 655 656 def create_default_functions(model_info): 657 """ 658 Autogenerate missing functions, such as Iqxy from Iq. 659 660 This only works for Iqxy when Iq is written in python. :func:`make_source` 661 performs a similar role for Iq written in C. 662 """ 663 if model_info['Iq'] is not None and model_info['Iqxy'] is None: 664 partable = model_info['parameters'] 665 if partable.type['1d'] != partable.type['2d']: 666 raise ValueError("Iqxy model is missing") 667 Iq = model_info['Iq'] 668 def Iqxy(qx, qy, **kw): 669 return Iq(np.sqrt(qx**2 + qy**2), **kw) 670 model_info['Iqxy'] = Iqxy 671 654 672 655 673 def make_model_info(kernel_module): … … 678 696 model variants (e.g., the list of cases) or is None if there are no 679 697 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 698 * *par_type* categorizes the model parameters. See 685 699 :func:`categorize_parameters` for details. 686 700 * *demo* contains the *{parameter: value}* map used in compare (and maybe … … 699 713 *model_info* blocks for the composition objects. This allows us to 700 714 build complete product and mixture models from just the info. 715 * *max_pd* is the max polydispersity dimension. This is constant and 716 should not be reset. You may be able to change it when the program 717 starts by setting *sasmodels.generate.MAX_PD*. 701 718 702 719 """ 703 720 # TODO: maybe turn model_info into a class ModelDefinition 704 parameters = COMMON_PARAMETERS + kernel_module.parameters 721 #print("make parameter table", kernel_module.parameters) 722 parameters = make_parameter_table(kernel_module.parameters) 705 723 filename = abspath(kernel_module.__file__) 706 724 kernel_id = splitext(basename(filename))[0] … … 731 749 functions = "ER VR form_volume Iq Iqxy shape sesans".split() 732 750 model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 751 create_default_functions(model_info) 752 # Precalculate the monodisperse parameters 753 # TODO: make this a lazy evaluator 754 # make_model_info is called for every model on sasview startup 755 model_info['mono_details'] = mono_details(model_info) 733 756 return model_info 734 757 -
sasmodels/kernelcl.py
r8e0d974 rfec69dd 75 75 76 76 77 # Pragmas for enable OpenCL features. Be sure to protect them so that they 78 # still compile even if OpenCL is not present. 79 _F16_PRAGMA = """\ 80 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 81 # pragma OPENCL EXTENSION cl_khr_fp16: enable 82 #endif 83 """ 84 85 _F64_PRAGMA = """\ 86 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 87 # pragma OPENCL EXTENSION cl_khr_fp64: enable 88 #endif 89 """ 90 91 77 92 ENV = None 78 93 def environment(): … … 142 157 raise RuntimeError("%s not supported for devices"%dtype) 143 158 144 source = generate.convert_type(source, dtype) 159 source_list = [generate.convert_type(source, dtype)] 160 161 if dtype == generate.F16: 162 source_list.insert(0, _F16_PRAGMA) 163 elif dtype == generate.F64: 164 source_list.insert(0, _F64_PRAGMA) 165 145 166 # Note: USE_SINCOS makes the intel cpu slower under opencl 146 167 if context.devices[0].type == cl.device_type.GPU: 147 source = "#define USE_SINCOS\n" + source168 source_list.insert(0, "#define USE_SINCOS\n") 148 169 options = (get_fast_inaccurate_build_options(context.devices[0]) 149 170 if fast else []) 171 source = "\n".join(source) 150 172 program = cl.Program(context, source).build(options=options) 151 173 return program … … 220 242 key = "%s-%s-%s"%(name, dtype, fast) 221 243 if key not in self.compiled: 222 #print("compiling",name)244 print("compiling",name) 223 245 dtype = np.dtype(dtype) 224 program = compile_model(self.get_context(dtype), source, dtype, fast) 246 program = compile_model(self.get_context(dtype), 247 str(source), dtype, fast) 225 248 self.compiled[key] = program 226 249 return self.compiled[key] … … 314 337 self.program = None 315 338 316 def __call__(self, q_vectors):339 def make_kernel(self, q_vectors): 317 340 if self.program is None: 318 341 compiler = environment().compile_program 319 self.program = compiler(self.info['name'], self.source, self.dtype,320 self. fast)342 self.program = compiler(self.info['name'], self.source, 343 self.dtype, self.fast) 321 344 is_2d = len(q_vectors) == 2 322 345 kernel_name = generate.kernel_name(self.info, is_2d) … … 365 388 # at this point, so instead using 32, which is good on the set of 366 389 # architectures tested so far. 367 self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 390 if self.is_2d: 391 # Note: 17 rather than 15 because results is 2 elements 392 # longer than input. 393 width = ((self.nq+17)//16)*16 394 self.q = np.empty((width, 2), dtype=dtype) 395 self.q[:self.nq, 0] = q_vectors[0] 396 self.q[:self.nq, 1] = q_vectors[1] 397 else: 398 # Note: 33 rather than 31 because results is 2 elements 399 # longer than input. 400 width = ((self.nq+33)//32)*32 401 self.q = np.empty(width, dtype=dtype) 402 self.q[:self.nq] = q_vectors[0] 403 self.global_size = [self.q.shape[0]] 368 404 context = env.get_context(self.dtype) 369 self.global_size = [self.q_vectors[0].size]370 405 #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 ] 406 self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 407 hostbuf=self.q) 375 408 376 409 def release(self): … … 378 411 Free the memory. 379 412 """ 380 for b in self.q_buffers:381 b.release()382 self.q_buffers = []413 if self.q is not None: 414 self.q.release() 415 self.q = None 383 416 384 417 def __del__(self): … … 406 439 """ 407 440 def __init__(self, kernel, model_info, q_vectors, dtype): 441 max_pd = self.info['max_pd'] 442 npars = len(model_info['parameters'])-2 408 443 q_input = GpuInput(q_vectors, dtype) 444 self.dtype = dtype 445 self.dim = '2d' if q_input.is_2d else '1d' 409 446 self.kernel = kernel 410 447 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] 448 self.pd_stop_index = 4*max_pd-1 449 # plus three for the normalization values 450 self.result = np.empty(q_input.nq+3, q_input.dtype) 415 451 416 452 # Inputs and outputs for each kernel call … … 418 454 env = environment() 419 455 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, 456 457 # details is int32 data, padded to an 8 integer boundary 458 size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 459 self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 423 460 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):461 self.q_input = q_input # allocated by GpuInput above 462 463 self._need_release = [ self.result_b, self.q_input ] 464 465 def __call__(self, details, weights, values, cutoff): 429 466 real = (np.float32 if self.q_input.dtype == generate.F32 430 467 else np.float64 if self.q_input.dtype == generate.F64 431 468 else np.float16 if self.q_input.dtype == generate.F16 432 469 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 470 assert details.dtype == np.int32 471 assert weights.dtype == real and values.dtype == real 472 473 context = self.queue.context 474 details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 475 hostbuf=details) 476 weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 477 hostbuf=weights) 478 values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 479 hostbuf=values) 480 481 start, stop = 0, self.details[self.pd_stop_index] 482 args = [ 483 np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 484 self.details_b, self.weights_b, self.values_b, 485 self.q_input.q_b, self.result_b, real(cutoff), 486 ] 460 487 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 488 cl.enqueue_copy(self.queue, self.result, self.result_b) 489 [v.release() for v in details_b, weights_b, values_b] 490 491 return self.result[:self.nq] 464 492 465 493 def release(self): -
sasmodels/kerneldll.py
r6ad0e87 r69aa451 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" … … 141 140 142 141 source = generate.convert_type(source, dtype) 143 source_files = generate.model_sources(model_info) + [model_info['filename']]142 newest = generate.timestamp(model_info) 144 143 dll = dll_path(model_info, dtype) 145 newest = max(os.path.getmtime(f) for f in source_files)146 144 if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 147 145 # Replace with a proper temp file … … 171 169 return DllModel(filename, model_info, dtype=dtype) 172 170 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 171 class DllModel(object): 178 172 """ … … 197 191 198 192 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 193 #print("dll", self.dllpath) 205 194 try: … … 212 201 else c_double if self.dtype == generate.F64 213 202 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 [] 203 204 # int, int, int, int*, double*, double*, double*, double*, double*, double 205 argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 216 206 self.Iq = self.dll[generate.kernel_name(self.info, False)] 217 self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d218 219 207 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() 208 self.Iq.argtypes = argtypes 209 self.Iqxy.argtypes = argtypes 223 210 224 211 def __getstate__(self): … … 229 216 self.dll = None 230 217 231 def __call__(self, q_vectors):218 def make_kernel(self, q_vectors): 232 219 q_input = PyInput(q_vectors, self.dtype) 233 220 if self.dll is None: self._load_dll() … … 272 259 """ 273 260 def __init__(self, kernel, model_info, q_input): 261 self.kernel = kernel 274 262 self.info = model_info 275 263 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): 264 self.dtype = q_input.dtype 265 self.dim = '2d' if q_input.is_2d else '1d' 266 self.result = np.empty(q_input.nq+3, q_input.dtype) 267 268 def __call__(self, details, weights, values, cutoff): 286 269 real = (np.float32 if self.q_input.dtype == generate.F32 287 270 else np.float64 if self.q_input.dtype == generate.F64 288 271 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) 272 assert details.dtype == np.int32 273 assert weights.dtype == real and values.dtype == real 274 275 max_pd = self.info['max_pd'] 276 start, stop = 0, details[4*max_pd-1] 277 print(details) 278 args = [ 279 self.q_input.nq, # nq 280 start, # pd_start 281 stop, # pd_stop pd_stride[MAX_PD] 282 details.ctypes.data, # problem 283 weights.ctypes.data, # weights 284 values.ctypes.data, #pars 285 self.q_input.q.ctypes.data, #q 286 self.result.ctypes.data, # results 287 real(cutoff), # cutoff 288 ] 303 289 self.kernel(*args) 304 305 return self.res 290 return self.result[:-3] 306 291 307 292 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/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/product.py
r8e0d974 r8e0d974 99 99 # a parameter map. 100 100 par_map = {} 101 p_info = p_kernel.info['par type']102 s_info = s_kernel.info['par type']101 p_info = p_kernel.info['par_type'] 102 s_info = s_kernel.info['par_type'] 103 103 vol_pars = set(p_info['volume']) 104 104 if dim == '2d': -
sasmodels/resolution.py
r8b935d1 r48fbd50 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)
Note: See TracChangeset
for help on using the changeset viewer.