Changeset 4937980 in sasmodels
- Timestamp:
- Apr 6, 2016 7:46:11 PM (9 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- cb3a824, 8bd7b77
- Parents:
- 541e7c9 (diff), e6408d0 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Files:
-
- 5 added
- 39 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/generate.py
rd2bb604 r4937980 15 15 16 16 *form_volume(p1, p2, ...)* returns the volume of the form with particular 17 dimension .17 dimension, or 1.0 if no volume normalization is required. 18 18 19 19 *ER(p1, p2, ...)* returns the effective radius of the form with -
sasmodels/models/flexible_cylinder.c
rd2bb604 r4937980 2 2 form_volume(length, kuhn_length, radius) 3 3 { 4 return 1. ;4 return 1.0; 5 5 } 6 6 -
sasmodels/models/lib/sas_JN.c
r541e7c9 r4937980 53 53 54 54 const double MACHEP = 1.11022302462515654042E-16; 55 double pkm2, pkm1, pk, xk, r, ans , xinv;55 double pkm2, pkm1, pk, xk, r, ans; 56 56 int k, sign; 57 57 … … 123 123 124 124 #else 125 xinv = 1.0/x;125 const double xinv = 1.0/x; 126 126 pkm1 = ans * xinv; 127 127 k = n-1; -
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 r21b116f 81 81 from bumps.names import Parameter 82 82 83 pars = {} 84 for p in model_info['parameters']: 83 pars = {} # floating point parameters 84 pd_types = {} # distribution names 85 for p in model_info['parameters'].call_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 name = p.name + '_pd_type' 98 pd_types[name] = kwargs.pop(name, 'gaussian') 102 99 103 100 if kwargs: # args not corresponding to parameters -
sasmodels/compare.py
rf247314 r21b116f 38 38 from . import core 39 39 from . import kerneldll 40 from . import product41 40 from .data import plot_theory, empty_data1D, empty_data2D 42 41 from .direct_model import DirectModel … … 306 305 Format the parameter list for printing. 307 306 """ 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 307 lines = [] 315 for p in model_info['parameters']:316 if exclude(p.name): continue308 parameters = model_info['parameters'] 309 for p in parameters.user_parameters(pars, is2d): 317 310 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'))311 value=pars.get(p.id, p.default), 312 pd=pars.get(p.id+"_pd", 0.), 313 n=int(pars.get(p.id+"_pd_n", 0)), 314 nsigma=pars.get(p.id+"_pd_nsgima", 3.), 315 type=pars.get(p.id+"_pd_type", 'gaussian')) 323 316 lines.append(_format_par(p.name, **fields)) 324 317 return "\n".join(lines) … … 454 447 """ 455 448 # initialize the code so time is more accurate 456 value = calculator(**suppress_pd(pars)) 449 if Nevals > 1: 450 value = calculator(**suppress_pd(pars)) 457 451 toc = tic() 458 452 for _ in range(max(Nevals, 1)): # make sure there is at least one eval … … 661 655 """ 662 656 # 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" 657 pars = {} 658 for p in model_info['parameters'].call_parameters: 659 parts = [('', p.default)] 660 if p.polydisperse: 661 parts.append(('_pd', 0.0)) 662 parts.append(('_pd_n', 0)) 663 parts.append(('_pd_nsigma', 3.0)) 664 parts.append(('_pd_type', "gaussian")) 665 for ext, val in parts: 666 if p.length > 1: 667 dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length)) 668 else: 669 pars[p.id+ext] = val 672 670 673 671 # Plug in values given in demo … … 774 772 775 773 if len(engines) == 0: 776 engines.extend(['single', ' sasview'])774 engines.extend(['single', 'double']) 777 775 elif len(engines) == 1: 778 if engines[0][0] != ' sasview':779 engines.append(' sasview')776 if engines[0][0] != 'double': 777 engines.append('double') 780 778 else: 781 779 engines.append('single') … … 879 877 model_info = opts['def'] 880 878 pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 879 # Initialize parameter ranges, fixing the 2D parameters for 1D data. 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'].user_parameters(is2d=False): 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/compare_many.py
rce346b6 ree8f734 20 20 21 21 from . import core 22 from . import generate23 22 from .compare import (MODELS, randomize_pars, suppress_pd, make_data, 24 23 make_engine, get_pars, columnize, … … 125 124 try: 126 125 result = fn(**pars) 127 except KeyboardInterrupt: 128 raise 129 except: 126 except Exception: 130 127 traceback.print_exc() 131 128 print("when comparing %s for %d"%(name, seed)) … … 246 243 base = sys.argv[5] 247 244 comp = sys.argv[6] if len(sys.argv) > 6 else "sasview" 248 except :245 except Exception: 249 246 traceback.print_exc() 250 247 print_usage() -
sasmodels/convert.json
rec45c4f rd2bb604 623 623 "SphereSLDModel", 624 624 { 625 "n": "n_shells", 625 626 "radius_core": "rad_core", 626 627 "sld_solvent": "sld_solv" -
sasmodels/core.py
r4d76711 ree8f734 24 24 from . import kernelcl 25 25 HAVE_OPENCL = True 26 except :26 except Exception: 27 27 HAVE_OPENCL = False 28 28 … … 64 64 """ 65 65 try: s + '' 66 except : return False66 except Exception: return False 67 67 return True 68 68 … … 170 170 171 171 172 def get_weights( model_info, pars, name):172 def get_weights(parameter, values): 173 173 """ 174 174 Generate the distribution for parameter *name* given the parameter values … … 178 178 from the *pars* dictionary for parameter value and parameter dispersion. 179 179 """ 180 relative = name in model_info['partype']['pd-rel'] 181 limits = model_info['limits'][name] 182 disperser = pars.get(name+'_pd_type', 'gaussian') 183 value = pars.get(name, model_info['defaults'][name]) 184 npts = pars.get(name+'_pd_n', 0) 185 width = pars.get(name+'_pd', 0.0) 186 nsigma = pars.get(name+'_pd_nsigma', 3.0) 180 value = float(values.get(parameter.name, parameter.default)) 181 relative = parameter.relative_pd 182 limits = parameter.limits 183 disperser = values.get(parameter.name+'_pd_type', 'gaussian') 184 npts = values.get(parameter.name+'_pd_n', 0) 185 width = values.get(parameter.name+'_pd', 0.0) 186 nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 187 if npts == 0 or width == 0: 188 return [value], [] 187 189 value, weight = weights.get_weights( 188 190 disperser, npts, width, nsigma, value, limits, relative) … … 198 200 """ 199 201 value, weight = zip(*pars) 202 weight = [w if w else [1.] for w in weight] 200 203 value = [v.flatten() for v in meshgrid(*value)] 201 204 weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) … … 216 219 *mono* is True if polydispersity should be set to none on all parameters. 217 220 """ 218 fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 219 for name in kernel.fixed_pars] 221 parameters = kernel.info['parameters'] 220 222 if mono: 221 pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 222 for name in kernel.pd_pars] 223 else: 224 pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 225 return kernel(fixed_pars, pd_pars, cutoff=cutoff) 223 active = lambda name: False 224 elif kernel.dim == '1d': 225 active = lambda name: name in parameters.pd_1d 226 elif kernel.dim == '2d': 227 active = lambda name: name in parameters.pd_2d 228 else: 229 active = lambda name: True 230 231 vw_pairs = [(get_weights(p, pars) if active(p.name) 232 else ([pars.get(p.name, p.default)], [])) 233 for p in parameters.call_parameters] 234 235 details, weights, values = build_details(kernel, vw_pairs) 236 return kernel(details, weights, values, cutoff) 237 238 def build_details(kernel, pairs): 239 values, weights = zip(*pairs) 240 if max([len(w) for w in weights]) > 1: 241 details = generate.poly_details(kernel.info, weights) 242 else: 243 details = kernel.info['mono_details'] 244 weights, values = [np.hstack(v) for v in (weights, values)] 245 weights = weights.astype(dtype=kernel.dtype) 246 values = values.astype(dtype=kernel.dtype) 247 return details, weights, values 248 226 249 227 250 def call_ER_VR(model_info, vol_pars): … … 256 279 return 1.0 257 280 else: 258 vol_pars = [get_weights(model_info, values, name) 259 for name in model_info['partype']['volume']] 281 vol_pars = [get_weights(parameter, values) 282 for parameter in model_info['parameters'].call_parameters 283 if parameter.type == 'volume'] 260 284 value, weight = dispersion_mesh(vol_pars) 261 285 individual_radii = ER(*value) 262 #print(values[0].shape, weights.shape, fv.shape)263 286 return np.sum(weight*individual_radii) / np.sum(weight) 264 287 … … 273 296 return 1.0 274 297 else: 275 vol_pars = [get_weights(model_info, values, name) 276 for name in model_info['partype']['volume']] 298 vol_pars = [get_weights(parameter, values) 299 for parameter in model_info['parameters'].call_parameters 300 if parameter.type == 'volume'] 277 301 value, weight = dispersion_mesh(vol_pars) 278 302 whole, part = VR(*value) -
sasmodels/data.py
rd6f5da6 ree8f734 358 358 try: 359 359 return fn(*args, **kw) 360 except KeyboardInterrupt: 361 raise 362 except: 360 except Exception: 363 361 traceback.print_exc() 364 362 -
sasmodels/direct_model.py
r4d76711 ree8f734 67 67 self.data_type = 'Iq' 68 68 69 partype = model.info['partype']70 71 69 if self.data_type == 'sesans': 72 70 … … 82 80 q_mono = sesans.make_all_q(data) 83 81 elif self.data_type == 'Iqxy': 84 if not partype['orientation'] and not partype['magnetic']:85 raise ValueError("not 2D without orientation or magnetic parameters")82 #if not model.info['parameters'].has_2d: 83 # raise ValueError("not 2D without orientation or magnetic parameters") 86 84 q = np.sqrt(data.qx_data**2 + data.qy_data**2) 87 85 qmin = getattr(data, 'qmin', 1e-16) … … 172 170 def _calc_theory(self, pars, cutoff=0.0): 173 171 if self._kernel is None: 174 self._kernel = self._model.make_kernel(self._kernel_inputs) # pylint: disable=attribute-dedata_type172 self._kernel = self._model.make_kernel(self._kernel_inputs) 175 173 self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs) 176 174 if self._kernel_mono_inputs else None) … … 243 241 try: 244 242 values = [float(v) for v in call.split(',')] 245 except :243 except Exception: 246 244 values = [] 247 245 if len(values) == 1: -
sasmodels/kernelcl.py
r4d76711 rd2bb604 48 48 harmless, albeit annoying. 49 49 """ 50 from __future__ import print_function 50 51 import os 51 52 import warnings … … 54 55 55 56 try: 57 raise NotImplementedError("OpenCL not yet implemented for new kernel template") 56 58 import pyopencl as cl 57 59 # Ask OpenCL for the default context so that we know that one exists … … 73 75 # of polydisperse parameters. 74 76 MAX_LOOPS = 2048 77 78 79 # Pragmas for enable OpenCL features. Be sure to protect them so that they 80 # still compile even if OpenCL is not present. 81 _F16_PRAGMA = """\ 82 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 83 # pragma OPENCL EXTENSION cl_khr_fp16: enable 84 #endif 85 """ 86 87 _F64_PRAGMA = """\ 88 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 89 # pragma OPENCL EXTENSION cl_khr_fp64: enable 90 #endif 91 """ 75 92 76 93 … … 142 159 raise RuntimeError("%s not supported for devices"%dtype) 143 160 144 source = generate.convert_type(source, dtype) 161 source_list = [generate.convert_type(source, dtype)] 162 163 if dtype == generate.F16: 164 source_list.insert(0, _F16_PRAGMA) 165 elif dtype == generate.F64: 166 source_list.insert(0, _F64_PRAGMA) 167 145 168 # Note: USE_SINCOS makes the intel cpu slower under opencl 146 169 if context.devices[0].type == cl.device_type.GPU: 147 source = "#define USE_SINCOS\n" + source170 source_list.insert(0, "#define USE_SINCOS\n") 148 171 options = (get_fast_inaccurate_build_options(context.devices[0]) 149 172 if fast else []) 173 source = "\n".join(source_list) 150 174 program = cl.Program(context, source).build(options=options) 151 175 return program … … 220 244 key = "%s-%s-%s"%(name, dtype, fast) 221 245 if key not in self.compiled: 222 #print("compiling",name)246 print("compiling",name) 223 247 dtype = np.dtype(dtype) 224 program = compile_model(self.get_context(dtype), source, dtype, fast) 248 program = compile_model(self.get_context(dtype), 249 str(source), dtype, fast) 225 250 self.compiled[key] = program 226 251 return self.compiled[key] … … 317 342 if self.program is None: 318 343 compiler = environment().compile_program 319 self.program = compiler(self.info['name'], self.source, self.dtype,320 self. fast)344 self.program = compiler(self.info['name'], self.source, 345 self.dtype, self.fast) 321 346 is_2d = len(q_vectors) == 2 322 347 kernel_name = generate.kernel_name(self.info, is_2d) … … 365 390 # at this point, so instead using 32, which is good on the set of 366 391 # architectures tested so far. 367 self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 392 if self.is_2d: 393 # Note: 17 rather than 15 because results is 2 elements 394 # longer than input. 395 width = ((self.nq+17)//16)*16 396 self.q = np.empty((width, 2), dtype=dtype) 397 self.q[:self.nq, 0] = q_vectors[0] 398 self.q[:self.nq, 1] = q_vectors[1] 399 else: 400 # Note: 33 rather than 31 because results is 2 elements 401 # longer than input. 402 width = ((self.nq+33)//32)*32 403 self.q = np.empty(width, dtype=dtype) 404 self.q[:self.nq] = q_vectors[0] 405 self.global_size = [self.q.shape[0]] 368 406 context = env.get_context(self.dtype) 369 self.global_size = [self.q_vectors[0].size]370 407 #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 ] 408 self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 409 hostbuf=self.q) 375 410 376 411 def release(self): … … 378 413 Free the memory. 379 414 """ 380 for b in self.q_buffers:381 b.release()382 self.q_buffers = []415 if self.q is not None: 416 self.q.release() 417 self.q = None 383 418 384 419 def __del__(self): … … 406 441 """ 407 442 def __init__(self, kernel, model_info, q_vectors, dtype): 443 max_pd = model_info['max_pd'] 444 npars = len(model_info['parameters'])-2 408 445 q_input = GpuInput(q_vectors, dtype) 446 self.dtype = dtype 447 self.dim = '2d' if q_input.is_2d else '1d' 409 448 self.kernel = kernel 410 449 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] 450 self.pd_stop_index = 4*max_pd-1 451 # plus three for the normalization values 452 self.result = np.empty(q_input.nq+3, q_input.dtype) 415 453 416 454 # Inputs and outputs for each kernel call … … 418 456 env = environment() 419 457 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, 458 459 # details is int32 data, padded to an 8 integer boundary 460 size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 461 self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 423 462 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):463 self.q_input = q_input # allocated by GpuInput above 464 465 self._need_release = [ self.result_b, self.q_input ] 466 467 def __call__(self, details, weights, values, cutoff): 429 468 real = (np.float32 if self.q_input.dtype == generate.F32 430 469 else np.float64 if self.q_input.dtype == generate.F64 431 470 else np.float16 if self.q_input.dtype == generate.F16 432 471 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 472 assert details.dtype == np.int32 473 assert weights.dtype == real and values.dtype == real 474 475 context = self.queue.context 476 details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 477 hostbuf=details) 478 weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 479 hostbuf=weights) 480 values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 481 hostbuf=values) 482 483 start, stop = 0, self.details[self.pd_stop_index] 484 args = [ 485 np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 486 self.details_b, self.weights_b, self.values_b, 487 self.q_input.q_b, self.result_b, real(cutoff), 488 ] 460 489 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 490 cl.enqueue_copy(self.queue, self.result, self.result_b) 491 [v.release() for v in details_b, weights_b, values_b] 492 493 return self.result[:self.nq] 464 494 465 495 def release(self): -
sasmodels/kerneldll.py
r4d76711 r1e2a1ba 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): … … 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
r4d76711 r1e2a1ba 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/list_pars.py
ra84a0ca r21b116f 25 25 for name in sorted(MODELS): 26 26 model_info = load_model_info(name) 27 for p in model_info['parameters'] :27 for p in model_info['parameters'].kernel_parameters: 28 28 partable.setdefault(p.name, []) 29 29 partable[p.name].append(name) -
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
r4d76711 rd2bb604 89 89 90 90 if is_py: # kernel implemented in python 91 continue # TODO: re-enable python tests 91 92 test_name = "Model: %s, Kernel: python"%model_name 92 93 test_method_name = "test_%s_python" % model_name -
sasmodels/models/cylinder.c
r26141cb re9b1663d 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/hardsphere.py
rec45c4f rd2bb604 149 149 """ 150 150 151 Iqxy = """152 // never called since no orientation or magnetic parameters.153 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);154 """155 156 151 # ER defaults to 0.0 157 152 # VR defaults to 1.0 -
sasmodels/models/hayter_msa.py
rec45c4f rd2bb604 87 87 return 1.0; 88 88 """ 89 Iqxy = """90 // never called since no orientation or magnetic parameters.91 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);92 """93 89 # ER defaults to 0.0 94 90 # VR defaults to 1.0 -
sasmodels/models/lamellar.py
rec45c4f rd2bb604 82 82 """ 83 83 84 Iqxy = """85 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);86 """87 88 84 # ER defaults to 0.0 89 85 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg.py
rec45c4f rd2bb604 101 101 """ 102 102 103 Iqxy = """104 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);105 """106 107 103 # ER defaults to 0.0 108 104 # VR defaults to 1.0 -
sasmodels/models/lamellar_hg_stack_caille.py
rec45c4f rd2bb604 120 120 """ 121 121 122 Iqxy = """123 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);124 """125 126 122 # ER defaults to 0.0 127 123 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_caille.py
rec45c4f rd2bb604 104 104 """ 105 105 106 Iqxy = """107 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);108 """109 110 106 # ER defaults to 0.0 111 107 # VR defaults to 1.0 -
sasmodels/models/lamellar_stack_paracrystal.py
rec45c4f rd2bb604 132 132 """ 133 133 134 Iqxy = """135 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);136 """137 138 134 # ER defaults to 0.0 139 135 # VR defaults to 1.0 -
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
rec45c4f rea05c87 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.c
r13ed84c rd2bb604 1 1 double Iq(double q, double case_num, 2 double Na, double Phia, double va, double a_sld, double ba, 3 double Nb, double Phib, double vb, double b_sld, double bb, 4 double Nc, double Phic, double vc, double c_sld, double bc, 5 double Nd, double Phid, double vd, double d_sld, double bd, 2 double N[], double Phi[], double v[], double L[], double b[], 6 3 double Kab, double Kac, double Kad, 7 4 double Kbc, double Kbd, double Kcd 8 5 ); 9 6 10 double Iqxy(double qx, double qy, double case_num,11 double Na, double Phia, double va, double a_sld, double ba,12 double Nb, double Phib, double vb, double b_sld, double bb,13 double Nc, double Phic, double vc, double c_sld, double bc,14 double Nd, double Phid, double vd, double d_sld, double bd,15 double Kab, double Kac, double Kad,16 double Kbc, double Kbd, double Kcd17 );18 19 double form_volume(void);20 21 double form_volume(void)22 {23 return 1.0;24 }25 26 7 double Iq(double q, double case_num, 27 double Na, double Phia, double va, double La, double ba, 28 double Nb, double Phib, double vb, double Lb, double bb, 29 double Nc, double Phic, double vc, double Lc, double bc, 30 double Nd, double Phid, double vd, double Ld, double bd, 8 double N[], double Phi[], double v[], double L[], double b[], 31 9 double Kab, double Kac, double Kad, 32 10 double Kbc, double Kbd, double Kcd … … 36 14 #if 0 // Sasview defaults 37 15 if (icase <= 1) { 38 N a=Nb=1000.0;39 Phi a=Phib=0.0000001;16 N[0]=N[1]=1000.0; 17 Phi[0]=Phi[1]=0.0000001; 40 18 Kab=Kac=Kad=Kbc=Kbd=-0.0004; 41 L a=Lb=1e-12;42 v a=vb=100.0;43 b a=bb=5.0;19 L[0]=L[1]=1e-12; 20 v[0]=v[1]=100.0; 21 b[0]=b[1]=5.0; 44 22 } else if (icase <= 4) { 45 Phi a=0.0000001;23 Phi[0]=0.0000001; 46 24 Kab=Kac=Kad=-0.0004; 47 L a=1e-12;48 v a=100.0;49 b a=5.0;25 L[0]=1e-12; 26 v[0]=100.0; 27 b[0]=5.0; 50 28 } 51 29 #else 52 30 if (icase <= 1) { 53 N a=Nb=0.0;54 Phi a=Phib=0.0;31 N[0]=N[1]=0.0; 32 Phi[0]=Phi[1]=0.0; 55 33 Kab=Kac=Kad=Kbc=Kbd=0.0; 56 L a=Lb=Ld;57 v a=vb=vd;58 b a=bb=0.0;34 L[0]=L[1]=L[3]; 35 v[0]=v[1]=v[3]; 36 b[0]=b[1]=0.0; 59 37 } else if (icase <= 4) { 60 N a= 0.0;61 Phi a=0.0;38 N[0] = 0.0; 39 Phi[0]=0.0; 62 40 Kab=Kac=Kad=0.0; 63 L a=Ld;64 v a=vd;65 b a=0.0;41 L[0]=L[3]; 42 v[0]=v[3]; 43 b[0]=0.0; 66 44 } 67 45 #endif 68 46 69 const double Xa = q*q*b a*ba*Na/6.0;70 const double Xb = q*q*b b*bb*Nb/6.0;71 const double Xc = q*q*b c*bc*Nc/6.0;72 const double Xd = q*q*b d*bd*Nd/6.0;47 const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 48 const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 49 const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 50 const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 73 51 74 52 // limit as Xa goes to 0 is 1 … … 98 76 #if 0 99 77 const double S0aa = icase<5 100 ? 1.0 : N a*Phia*va*Paa;78 ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 101 79 const double S0bb = icase<2 102 ? 1.0 : N b*Phib*vb*Pbb;103 const double S0cc = N c*Phic*vc*Pcc;104 const double S0dd = N d*Phid*vd*Pdd;80 ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 81 const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 82 const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 105 83 const double S0ab = icase<8 106 ? 0.0 : sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;84 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 107 85 const double S0ac = icase<9 108 ? 0.0 : sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);86 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 109 87 const double S0ad = icase<9 110 ? 0.0 : sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);88 ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 111 89 const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 112 ? 0.0 : sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;90 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 113 91 const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 114 ? 0.0 : sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);92 ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 115 93 const double S0cd = (icase==0 || icase==2 || icase==5) 116 ? 0.0 : sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;94 ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 117 95 #else // sasview equivalent 118 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N c,Phic,vc,Pcc);119 double S0aa = N a*Phia*va*Paa;120 double S0bb = N b*Phib*vb*Pbb;121 double S0cc = N c*Phic*vc*Pcc;122 double S0dd = N d*Phid*vd*Pdd;123 double S0ab = sqrt(N a*va*Phia*Nb*vb*Phib)*Pa*Pb;124 double S0ac = sqrt(N a*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);125 double S0ad = sqrt(N a*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);126 double S0bc = sqrt(N b*vb*Phib*Nc*vc*Phic)*Pb*Pc;127 double S0bd = sqrt(N b*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);128 double S0cd = sqrt(N c*vc*Phic*Nd*vd*Phid)*Pc*Pd;96 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 97 double S0aa = N[0]*Phi[0]*v[0]*Paa; 98 double S0bb = N[1]*Phi[1]*v[1]*Pbb; 99 double S0cc = N[2]*Phi[2]*v[2]*Pcc; 100 double S0dd = N[3]*Phi[3]*v[3]*Pdd; 101 double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 102 double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 103 double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 104 double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 105 double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 106 double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 129 107 switch(icase){ 130 108 case 0: … … 311 289 // Note: 1e-13 to convert from fm to cm for scattering length 312 290 const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 313 const double Lad = icase<5 ? 0.0 : (L a/va - Ld/vd)*sqrt_Nav;314 const double Lbd = icase<2 ? 0.0 : (L b/vb - Ld/vd)*sqrt_Nav;315 const double Lcd = (L c/vc - Ld/vd)*sqrt_Nav;291 const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 292 const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 293 const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 316 294 317 295 const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 … … 321 299 322 300 } 323 324 double Iqxy(double qx, double qy,325 double case_num,326 double Na, double Phia, double va, double a_sld, double ba,327 double Nb, double Phib, double vb, double b_sld, double bb,328 double Nc, double Phic, double vc, double c_sld, double bc,329 double Nd, double Phid, double vd, double d_sld, double bd,330 double Kab, double Kac, double Kad,331 double Kbc, double Kbd, double Kcd332 )333 {334 double q = sqrt(qx*qx + qy*qy);335 return Iq(q,336 case_num,337 Na, Phia, va, a_sld, ba,338 Nb, Phib, vb, b_sld, bb,339 Nc, Phic, vc, c_sld, bc,340 Nd, Phid, vd, d_sld, bd,341 Kab, Kac, Kad,342 Kbc, Kbd, Kcd);343 } -
sasmodels/models/rpa.py
rec45c4f rea05c87 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
rec45c4f rd2bb604 170 170 # pylint: disable=bad-whitespace, line-too-long 171 171 # ["name", "units", default, [lower, upper], "type", "description"], 172 parameters = [["n _shells","", 1, [0, 9], "", "number of shells"],172 parameters = [["n", "", 1, [0, 9], "", "number of shells"], 173 173 ["radius_core", "Ang", 50.0, [0, inf], "", "intern layer thickness"], 174 174 ["sld_core", "1e-6/Ang^2", 2.07, [-inf, inf], "", "sld function flat"], … … 192 192 193 193 demo = dict( 194 n _shells=4,194 n=4, 195 195 scale=1.0, 196 196 solvent_sld=1.0, -
sasmodels/models/squarewell.py
rec45c4f rd2bb604 123 123 """ 124 124 125 Iqxy = """126 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);127 """128 129 125 # ER defaults to 0.0 130 126 # VR defaults to 1.0 -
sasmodels/models/stickyhardsphere.py
rec45c4f rd2bb604 171 171 """ 172 172 173 Iqxy = """174 return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);175 """176 177 173 # ER defaults to 0.0 178 174 # VR defaults to 1.0 -
sasmodels/product.py
rf247314 rea05c87 14 14 15 15 from .core import call_ER_VR 16 from .generate import process_parameters17 16 18 17 SCALE=0 … … 58 57 # Iq, Iqxy, form_volume, ER, VR and sesans 59 58 model_info['composition'] = ('product', [p_info, s_info]) 60 process_parameters(model_info)61 59 return model_info 62 60 … … 93 91 # a parameter map. 94 92 par_map = {} 95 p_info = p_kernel.info['par type']96 s_info = s_kernel.info['par type']93 p_info = p_kernel.info['par_type'] 94 s_info = s_kernel.info['par_type'] 97 95 vol_pars = set(p_info['volume']) 98 96 if dim == '2d': -
sasmodels/resolution.py
r4d76711 rd2bb604 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()) - keys 507 raise ValueError("bad parameters: [%s] not in [%s]"% 508 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 504 par_set = set([p.name for p in form.info['parameters'].call_parameters]) 505 if any(k not in par_set for k in pars.keys()): 506 extra = set(pars.keys()) - par_set 507 raise ValueError("bad parameters: [%s] not in [%s]" 508 % (", ".join(sorted(extra)), 509 ", ".join(sorted(pars.keys())))) 509 510 510 511 if np.isscalar(width): … … 556 557 from scipy.integrate import romberg 557 558 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()) - keys 561 raise ValueError("bad parameters: [%s] not in [%s]"% 562 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 559 par_set = set([p.name for p in form.info['parameters'].call_parameters]) 560 if any(k not in par_set for k in pars.keys()): 561 extra = set(pars.keys()) - par_set 562 raise ValueError("bad parameters: [%s] not in [%s]" 563 % (", ".join(sorted(extra)), 564 ", ".join(sorted(pars.keys())))) 563 565 564 566 _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) -
sasmodels/sasview_model.py
r4d76711 ree8f734 26 26 from . import custom 27 27 from . import generate 28 from . import weights 28 29 29 30 def load_standard_models(): … … 38 39 try: 39 40 models.append(_make_standard_model(name)) 40 except :41 except Exception: 41 42 logging.error(traceback.format_exc()) 42 43 return models … … 84 85 self._model = None 85 86 model_info = self._model_info 87 parameters = model_info['parameters'] 86 88 87 89 self.name = model_info['name'] 88 90 self.description = model_info['description'] 89 91 self.category = None 90 self.multiplicity_info = None 91 self.is_multifunc = False 92 #self.is_multifunc = False 93 for p in parameters.kernel_parameters: 94 if p.is_control: 95 profile_axes = model_info['profile_axes'] 96 self.multiplicity_info = [ 97 p.limits[1], p.name, p.choices, profile_axes[0] 98 ] 99 break 100 else: 101 self.multiplicity_info = [] 92 102 93 103 ## interpret the parameters … … 96 106 self.params = collections.OrderedDict() 97 107 self.dispersion = dict() 98 partype = model_info['partype'] 99 100 for p in model_info['parameters']: 108 109 self.orientation_params = [] 110 self.magnetic_params = [] 111 self.fixed = [] 112 for p in parameters.user_parameters(): 101 113 self.params[p.name] = p.default 102 114 self.details[p.name] = [p.units] + p.limits 103 104 for name in partype['pd-2d']: 105 self.dispersion[name] = { 106 'width': 0, 107 'npts': 35, 108 'nsigmas': 3, 109 'type': 'gaussian', 110 } 111 112 self.orientation_params = ( 113 partype['orientation'] 114 + [n + '.width' for n in partype['orientation']] 115 + partype['magnetic']) 116 self.magnetic_params = partype['magnetic'] 117 self.fixed = [n + '.width' for n in partype['pd-2d']] 115 if p.polydisperse: 116 self.dispersion[p.name] = { 117 'width': 0, 118 'npts': 35, 119 'nsigmas': 3, 120 'type': 'gaussian', 121 } 122 if p.type == 'orientation': 123 self.orientation_params.append(p.name) 124 self.orientation_params.append(p.name+".width") 125 self.fixed.append(p.name+".width") 126 if p.type == 'magnetic': 127 self.orientation_params.append(p.name) 128 self.magnetic_params.append(p.name) 129 self.fixed.append(p.name+".width") 130 118 131 self.non_fittable = [] 119 132 … … 234 247 """ 235 248 # TODO: fix test so that parameter order doesn't matter 236 ret = ['%s.%s' % (d.lower(), p) 237 for d in self._model_info['partype']['pd-2d'] 238 for p in ('npts', 'nsigmas', 'width')] 249 ret = ['%s.%s' % (p.name.lower(), ext) 250 for p in self._model_info['parameters'].user_parameters() 251 for ext in ('npts', 'nsigmas', 'width') 252 if p.polydisperse] 239 253 #print(ret) 240 254 return ret … … 309 323 # Check whether we have a list of ndarrays [qx,qy] 310 324 qx, qy = qdist 311 partype = self._model_info['partype'] 312 if not partype['orientation'] and not partype['magnetic']: 325 if not self._model_info['parameters'].has_2d: 313 326 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 314 327 else: … … 335 348 self._model = core.build_model(self._model_info) 336 349 q_vectors = [np.asarray(q) for q in args] 337 fn = self._model.make_kernel(q_vectors) 338 pars = [self.params[v] for v in fn.fixed_pars] 339 pd_pars = [self._get_weights(p) for p in fn.pd_pars] 340 result = fn(pars, pd_pars, self.cutoff) 341 fn.q_input.release() 342 fn.release() 350 kernel = self._model.make_kernel(q_vectors) 351 pairs = [self._get_weights(p) 352 for p in self._model_info['parameters'].call_parameters] 353 details, weights, values = core.build_details(kernel, pairs) 354 result = kernel(details, weights, values, cutoff=self.cutoff) 355 kernel.q_input.release() 356 kernel.release() 343 357 return result 344 358 … … 415 429 Return dispersion weights for parameter 416 430 """ 417 from . import weights 418 relative = self._model_info['partype']['pd-rel'] 419 limits = self._model_info['limits'] 420 dis = self.dispersion[par] 421 value, weight = weights.get_weights( 422 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 423 self.params[par], limits[par], par in relative) 424 return value, weight / np.sum(weight) 425 431 if par.polydisperse: 432 dis = self.dispersion[par.name] 433 value, weight = weights.get_weights( 434 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 435 self.params[par.name], par.limits, par.relative_pd) 436 return value, weight / np.sum(weight) 437 else: 438 return [self.params[par.name]], [] 426 439 427 440 def test_model():
Note: See TracChangeset
for help on using the changeset viewer.