Changeset 3543141 – SasView

sasmodels/models/broad_peak.py

rec45c4f	r65279d8
101	101	"""
102	102	return Iq(sqrt(qx 2 + qy 2), *args)
103
104	103	Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
105	104

sasmodels/models/poly_gauss_coil.py

-                      rec45c4f
+                      r65279d8
 """
+from numpy import inf, sqrt, exp, power
+import numpy as np
+from numpy import inf, exp, power, sqrt
 name =  "poly_gauss_coil"
 …
     # pylint: disable = missing-docstring
     u = polydispersity - 1.0
+    z = ((q * radius_gyration) * (q * radius_gyration)) / (1.0 + 2.0 * u)
+    if (q == 0).any():
+        inten = i_zero
+    z = (q*radius_gyration)**2 / (1.0 + 2.0*u)
+    # need to trap the case of the polydispersity being 1 (ie, monodispersity!)
+    if polydispersity == 1.0:
+        inten = i_zero * 2.0 * (exp(-z) + z - 1.0)
     else:
+    # need to trap the case of the polydispersity being 1 (ie, monodispersity!)
+        if polydispersity == 1:
+            inten = i_zero * 2.0 * (exp(-z) + z - 1.0 ) / (z * z)
+        else:
+            minusoneonu = -1.0 / u
+            inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z))
+        inten = i_zero * 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u)
+    index = q != 0.
+    inten[~index] = i_zero
+    inten[index] /= z[index]**2
     return inten
 #Iq.vectorized =  True  # Iq accepts an array of q values
+Iq.vectorized =  True  # Iq accepts an array of q values
 def Iqxy(qx, qy, *args):
     # pylint: disable = missing-docstring
     return Iq(sqrt(qx ** 2 + qy ** 2), *args)
 #Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
+Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values
 demo =  dict(scale = 1.0,

doc/developer/index.rst

-                      r56fc97a
+                      rb85be2d
 ***************************
+.. toctree::
+   :numbered: 4
+   :maxdepth: 4
+   calculator.rst

sasmodels/bumps_model.py

-                      rea75043
+                      r21b116f
     from bumps.names import Parameter
+    pars = {}
+    for p in model_info['parameters']:
+    pars = {}     # floating point parameters
+    pd_types = {} # distribution names
+    for p in model_info['parameters'].call_parameters:
         value = kwargs.pop(p.name, p.default)
         pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits)
+    for name in model_info['partype']['pd-2d']:
+        for xpart, xdefault, xlimits in [
+                ('_pd', 0., pars[name].limits),
+                ('_pd_n', 35., (0, 1000)),
+                ('_pd_nsigma', 3., (0, 10)),
+            ]:
+            xname = name + xpart
+            xvalue = kwargs.pop(xname, xdefault)
+            pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits)
+    pd_types = {}
+    for name in model_info['partype']['pd-2d']:
+        xname = name + '_pd_type'
+        xvalue = kwargs.pop(xname, 'gaussian')
+        pd_types[xname] = xvalue
+        if p.polydisperse:
+            for part, default, limits in [
+                    ('_pd', 0., pars[p.name].limits),
+                    ('_pd_n', 35., (0, 1000)),
+                    ('_pd_nsigma', 3., (0, 10)),
+                ]:
+                name = p.name + part
+                value = kwargs.pop(name, default)
+                pars[name] = Parameter.default(value, name=name, limits=limits)
+            name = p.name + '_pd_type'
+            pd_types[name] = kwargs.pop(name, 'gaussian')
     if kwargs:  # args not corresponding to parameters

sasmodels/compare.py

-                      rf247314
+                      r8bd7b77
 from . import core
 from . import kerneldll
-from . import product
 from .data import plot_theory, empty_data1D, empty_data2D
 from .direct_model import DirectModel
 …
     Choose a parameter range based on parameter name and initial value.
     """
+    # process the polydispersity options
     if p.endswith('_pd_n'):
         return [0, 100]
 …
         else:
             return [-180, 180]
+    elif p.endswith('_pd'):
+        return [0, 1]
     elif 'sld' in p:
         return [-0.5, 10]
-    elif p.endswith('_pd'):
-        return [0, 1]
     elif p == 'background':
         return [0, 10]
     elif p == 'scale':
         return [0, 1e3]
-    elif p == 'case_num':
-        # RPA hack
-        return [0, 10]
     elif v < 0:
-        # Kxy parameters in rpa model can be negative
         return [2*v, -2*v]
     else:
 …
 def _randomize_one(p, v):
+def _randomize_one(model_info, p, v):
     """
     Randomize a single parameter.
 …
     if any(p.endswith(s) for s in ('_pd_n', '_pd_nsigma', '_pd_type')):
         return v
+    # Find the parameter definition
+    for par in model_info['parameters'].call_parameters:
+        if par.name == p:
+            break
     else:
+        return np.random.uniform(*parameter_range(p, v))
+def randomize_pars(pars, seed=None):
+        raise ValueError("unknown parameter %r"%p)
+    if len(par.limits) > 2:  # choice list
+        return np.random.randint(len(par.limits))
+    limits = par.limits
+    if not np.isfinite(limits).all():
+        low, high = parameter_range(p, v)
+        limits = (max(limits[0], low), min(limits[1], high))
+    return np.random.uniform(*limits)
+def randomize_pars(model_info, pars, seed=None):
     """
     Generate random values for all of the parameters.
 …
     with push_seed(seed):
         # Note: the sort guarantees order `of calls to random number generator
         pars = dict((p, _randomize_one(p, v))
+        pars = dict((p, _randomize_one(model_info, p, v))
                     for p, v in sorted(pars.items()))
     return pars
 …
         # Make sure phi sums to 1.0
         if pars['case_num'] < 2:
             pars['Phia'] = 0.
             pars['Phib'] = 0.
+            pars['Phi1'] = 0.
+            pars['Phi2'] = 0.
         elif pars['case_num'] < 5:
             pars['Phia'] = 0.
         total = sum(pars['Phi'+c] for c in 'abcd')
         for c in 'abcd':
+            pars['Phi1'] = 0.
+        total = sum(pars['Phi'+c] for c in '1234')
+        for c in '1234':
             pars['Phi'+c] /= total
 …
     Format the parameter list for printing.
     """
-    if is2d:
-        exclude = lambda n: False
-    else:
-        partype = model_info['partype']
-        par1d = set(partype['fixed-1d']+partype['pd-1d'])
-        exclude = lambda n: n not in par1d
     lines = []
     for p in model_info['parameters']:
         if exclude(p.name): continue
+    parameters = model_info['parameters']
+    for p in parameters.user_parameters(pars, is2d):
         fields = dict(
             value=pars.get(p.name, p.default),
             pd=pars.get(p.name+"_pd", 0.),
             n=int(pars.get(p.name+"_pd_n", 0)),
             nsigma=pars.get(p.name+"_pd_nsgima", 3.),
             type=pars.get(p.name+"_pd_type", 'gaussian'))
+            value=pars.get(p.id, p.default),
+            pd=pars.get(p.id+"_pd", 0.),
+            n=int(pars.get(p.id+"_pd_n", 0)),
+            nsigma=pars.get(p.id+"_pd_nsgima", 3.),
+            type=pars.get(p.id+"_pd_type", 'gaussian'))
         lines.append(_format_par(p.name, **fields))
     return "\n".join(lines)
 …
     """
     # initialize the code so time is more accurate
+    value = calculator(**suppress_pd(pars))
+    if Nevals > 1:
+        value = calculator(**suppress_pd(pars))
     toc = tic()
     for _ in range(max(Nevals, 1)):  # make sure there is at least one eval
 …
     """
     # Get the default values for the parameters
+    pars = dict((p.name, p.default) for p in model_info['parameters'])
+    # Fill in default values for the polydispersity parameters
+    for p in model_info['parameters']:
+        if p.type in ('volume', 'orientation'):
+            pars[p.name+'_pd'] = 0.0
+            pars[p.name+'_pd_n'] = 0
+            pars[p.name+'_pd_nsigma'] = 3.0
+            pars[p.name+'_pd_type'] = "gaussian"
+    pars = {}
+    for p in model_info['parameters'].call_parameters:
+        parts = [('', p.default)]
+        if p.polydisperse:
+            parts.append(('_pd', 0.0))
+            parts.append(('_pd_n', 0))
+            parts.append(('_pd_nsigma', 3.0))
+            parts.append(('_pd_type', "gaussian"))
+        for ext, val in parts:
+            if p.length > 1:
+                dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length))
+            else:
+                pars[p.id+ext] = val
     # Plug in values given in demo
 …
     if len(engines) == 0:
         engines.extend(['single', 'sasview'])
+        engines.extend(['single', 'double'])
     elif len(engines) == 1:
         if engines[0][0] != 'sasview':
             engines.append('sasview')
+        if engines[0][0] != 'double':
+            engines.append('double')
         else:
             engines.append('single')
 …
     #pars.update(set_pars)  # set value before random to control range
     if opts['seed'] > -1:
         pars = randomize_pars(pars, seed=opts['seed'])
+        pars = randomize_pars(model_info, pars, seed=opts['seed'])
         print("Randomize using -random=%i"%opts['seed'])
     if opts['mono']:
 …
         model_info = opts['def']
         pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars'])
+        # Initialize parameter ranges, fixing the 2D parameters for 1D data.
         if not opts['is2d']:
+            active = [base + ext
+                      for base in model_info['partype']['pd-1d']
+                      for ext in ['', '_pd', '_pd_n', '_pd_nsigma']]
+            active.extend(model_info['partype']['fixed-1d'])
+            for k in active:
+                v = pars[k]
+                v.range(*parameter_range(k, v.value))
+            for p in model_info['parameters'].user_parameters(is2d=False):
+                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
+                    k = p.name+ext
+                    v = pars.get(k, None)
+                    if v is not None:
+                        v.range(*parameter_range(k, v.value))
         else:
             for k, v in pars.items():

sasmodels/compare_many.py

-                      rce346b6
+                      r8bd7b77
 from . import core
-from . import generate
 from .compare import (MODELS, randomize_pars, suppress_pd, make_data,
                       make_engine, get_pars, columnize,
 …
     if is_2d:
         if not model_info['has_2d']:
+        if not model_info['parameters'].has_2d:
             print(',"1-D only"')
             return
 …
         try:
             result = fn(**pars)
+        except KeyboardInterrupt:
+            raise
+        except:
+        except Exception:
             traceback.print_exc()
             print("when comparing %s for %d"%(name, seed))
 …
         base = sys.argv[5]
         comp = sys.argv[6] if len(sys.argv) > 6 else "sasview"
     except:
+    except Exception:
         traceback.print_exc()
         print_usage()

sasmodels/convert.json

-                      rec45c4f
+                      r6e7ff6d
     "RPAModel",
+    {
+      "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd",
+      "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld",
+      "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd",
+      "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd",
+      "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid",
       "case_num": "lcase_n"
+    }
 …
+    }
   ],
+  "_spherepy": [
+    "SphereModel",
+    {
+      "sld": "sldSph",
+      "radius": "radius",
+      "sld_solvent": "sldSolv"
+    }
+  ],
   "spherical_sld": [
     "SphereSLDModel",
+    {
+      "n": "n_shells",
       "radius_core": "rad_core",
       "sld_solvent": "sld_solv"

sasmodels/convert.py

rf247314	r8bd7b77
206	206	elif name == 'rpa':
207	207	# convert scattering lengths from femtometers to centimeters
208		for p in "L~~a", "Lb", "Lc", "Ld~~":
	208	for p in "L1", "L2", "L3", "L4":
209	209	if p in oldpars: oldpars[p] *= 1e-13
210	210	elif name == 'core_shell_parallelepiped':

sasmodels/core.py

-                      r4d76711
+                      ree8f734
     from . import kernelcl
     HAVE_OPENCL = True
 except:
+except Exception:
     HAVE_OPENCL = False
 …
     """
     try: s + ''
     except: return False
+    except Exception: return False
     return True
 …
 def get_weights(model_info, pars, name):
+def get_weights(parameter, values):
     """
     Generate the distribution for parameter *name* given the parameter values
 …
     from the *pars* dictionary for parameter value and parameter dispersion.
     """
+    relative = name in model_info['partype']['pd-rel']
+    limits = model_info['limits'][name]
+    disperser = pars.get(name+'_pd_type', 'gaussian')
+    value = pars.get(name, model_info['defaults'][name])
+    npts = pars.get(name+'_pd_n', 0)
+    width = pars.get(name+'_pd', 0.0)
+    nsigma = pars.get(name+'_pd_nsigma', 3.0)
+    value = float(values.get(parameter.name, parameter.default))
+    relative = parameter.relative_pd
+    limits = parameter.limits
+    disperser = values.get(parameter.name+'_pd_type', 'gaussian')
+    npts = values.get(parameter.name+'_pd_n', 0)
+    width = values.get(parameter.name+'_pd', 0.0)
+    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0)
+    if npts == 0 or width == 0:
+        return [value], []
     value, weight = weights.get_weights(
         disperser, npts, width, nsigma, value, limits, relative)
 …
     """
     value, weight = zip(*pars)
+    weight = [w if w else [1.] for w in weight]
     value = [v.flatten() for v in meshgrid(*value)]
     weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
 …
     *mono* is True if polydispersity should be set to none on all parameters.
     """
+    fixed_pars = [pars.get(name, kernel.info['defaults'][name])
+                  for name in kernel.fixed_pars]
+    parameters = kernel.info['parameters']
     if mono:
+        pd_pars = [( np.array([pars[name]]), np.array([1.0]) )
+                   for name in kernel.pd_pars]
+    else:
+        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars]
+    return kernel(fixed_pars, pd_pars, cutoff=cutoff)
+        active = lambda name: False
+    elif kernel.dim == '1d':
+        active = lambda name: name in parameters.pd_1d
+    elif kernel.dim == '2d':
+        active = lambda name: name in parameters.pd_2d
+    else:
+        active = lambda name: True
+    vw_pairs = [(get_weights(p, pars) if active(p.name)
+                 else ([pars.get(p.name, p.default)], []))
+                for p in parameters.call_parameters]
+    details, weights, values = build_details(kernel, vw_pairs)
+    return kernel(details, weights, values, cutoff)
+def build_details(kernel, pairs):
+    values, weights = zip(*pairs)
+    if max([len(w) for w in weights]) > 1:
+        details = generate.poly_details(kernel.info, weights)
+    else:
+        details = kernel.info['mono_details']
+    weights, values = [np.hstack(v) for v in (weights, values)]
+    weights = weights.astype(dtype=kernel.dtype)
+    values = values.astype(dtype=kernel.dtype)
+    return details, weights, values
 def call_ER_VR(model_info, vol_pars):
 …
         return 1.0
     else:
+        vol_pars = [get_weights(model_info, values, name)
+                    for name in model_info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters'].call_parameters
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         individual_radii = ER(*value)
-        #print(values[0].shape, weights.shape, fv.shape)
         return np.sum(weight*individual_radii) / np.sum(weight)
 …
         return 1.0
     else:
+        vol_pars = [get_weights(model_info, values, name)
+                    for name in model_info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters'].call_parameters
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         whole, part = VR(*value)

sasmodels/data.py

-                      rd6f5da6
+                      ree8f734
         try:
             return fn(*args, **kw)
+        except KeyboardInterrupt:
+            raise
+        except:
+        except Exception:
             traceback.print_exc()

sasmodels/direct_model.py

-                      r4d76711
+                      ree8f734
             self.data_type = 'Iq'
-        partype = model.info['partype']
         if self.data_type == 'sesans':
 …
             q_mono = sesans.make_all_q(data)
         elif self.data_type == 'Iqxy':
             if not partype['orientation'] and not partype['magnetic']:
                 raise ValueError("not 2D without orientation or magnetic parameters")
+            #if not model.info['parameters'].has_2d:
+            #    raise ValueError("not 2D without orientation or magnetic parameters")
             q = np.sqrt(data.qx_data**2 + data.qy_data**2)
             qmin = getattr(data, 'qmin', 1e-16)
 …
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
             self._kernel = self._model.make_kernel(self._kernel_inputs)  # pylint: disable=attribute-dedata_type
+            self._kernel = self._model.make_kernel(self._kernel_inputs)
             self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs)
                                  if self._kernel_mono_inputs else None)
 …
         try:
             values = [float(v) for v in call.split(',')]
         except:
+        except Exception:
             values = []
         if len(values) == 1:

sasmodels/generate.py

-                      re6408d0
+                      r6e7ff6d
     *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms.
+    #define INVALID(v) (expr)  returns False if v.parameter is invalid
+    for some parameter or other (e.g., v.bell_radius < v.radius).  If
+    necessary, the expression can call a function.
 These functions are defined in a kernel module .py script and an associated
 …
 The constructor code will generate the necessary vectors for computing
 them with the desired polydispersity.
-The available kernel parameters are defined as a list, with each parameter
-defined as a sublist with the following elements:
-    *name* is the name that will be used in the call to the kernel
-    function and the name that will be displayed to the user.  Names
-    should be lower case, with words separated by underscore.  If
-    acronyms are used, the whole acronym should be upper case.
-    *units* should be one of *degrees* for angles, *Ang* for lengths,
-    *1e-6/Ang^2* for SLDs.
-    *default value* will be the initial value for  the model when it
-    is selected, or when an initial value is not otherwise specified.
-    *limits = [lb, ub]* are the hard limits on the parameter value, used to
-    limit the polydispersity density function.  In the fit, the parameter limits
-    given to the fit are the limits  on the central value of the parameter.
-    If there is polydispersity, it will evaluate parameter values outside
-    the fit limits, but not outside the hard limits specified in the model.
-    If there are no limits, use +/-inf imported from numpy.
-    *type* indicates how the parameter will be used.  "volume" parameters
-    will be used in all functions.  "orientation" parameters will be used
-    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
-    *Imagnetic* only.  If *type* is the empty string, the parameter will
-    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters
-    can automatically be promoted to magnetic parameters, each of which
-    will have a magnitude and a direction, which may be different from
-    other sld parameters.
-    *description* is a short description of the parameter.  This will
-    be displayed in the parameter table and used as a tool tip for the
-    parameter value in the user interface.
 The kernel module must set variables defining the kernel meta data:
 …
 An *model_info* dictionary is constructed from the kernel meta data and
 returned to the caller.
-The model evaluator, function call sequence consists of q inputs and the return vector,
-followed by the loop value/weight vector, followed by the values for
-the non-polydisperse parameters, followed by the lengths of the
-polydispersity loops.  To construct the call for 1D models, the
-categories *fixed-1d* and *pd-1d* list the names of the parameters
-of the non-polydisperse and the polydisperse parameters respectively.
-Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models.
-The *pd-rel* category is a set of those parameters which give
-polydispersitiy as a portion of the value (so a 10% length dispersity
-would use a polydispersity value of 0.1) rather than absolute
-dispersity such as an angle plus or minus 15 degrees.
-The *volume* category lists the volume parameters in order for calls
-to volume within the kernel (used for volume normalization) and for
-calls to ER and VR for effective radius and volume ratio respectively.
-The *orientation* and *magnetic* categories list the orientation and
-magnetic parameters.  These are used by the sasview interface.  The
-blank category is for parameters such as scale which don't have any
-other marking.
 The doc string at the start of the kernel module will be used to
 …
 from __future__ import print_function
-#TODO: identify model files which have changed since loading and reload them.
 #TODO: determine which functions are useful outside of generate
 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
-import sys
 from os.path import abspath, dirname, join as joinpath, exists, basename, \
     splitext
+    splitext, getmtime
 import re
 import string
 import warnings
-from collections import namedtuple
 import numpy as np
+from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo
 from .custom import load_custom_kernel_module
+PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description']
+Parameter = namedtuple('Parameter', PARAMETER_FIELDS)
+C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
+TEMPLATE_ROOT = dirname(__file__)
 F16 = np.dtype('float16')
 …
 except TypeError:
     F128 = None
-# Scale and background, which are parameters common to every form factor
-COMMON_PARAMETERS = [
-    ["scale", "", 1, [0, np.inf], "", "Source intensity"],
-    ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"],
+    ]
 # Conversion from units defined in the parameter table for each model
 …
     raise ValueError("%r not found in %s" % (filename, search_path))
 def model_sources(model_info):
     """
 …
     return [_search(search_path, f) for f in model_info['source']]
+# Pragmas for enable OpenCL features.  Be sure to protect them so that they
+# still compile even if OpenCL is not present.
+_F16_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16)
+#  pragma OPENCL EXTENSION cl_khr_fp16: enable
+#endif
+"""
+_F64_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64)
+#  pragma OPENCL EXTENSION cl_khr_fp64: enable
+#endif
+"""
+def timestamp(model_info):
+    """
+    Return a timestamp for the model corresponding to the most recently
+    changed file or dependency.
+    """
+    source_files = (model_sources(model_info)
+                    + model_templates()
+                    + [model_info['filename']])
+    newest = max(getmtime(f) for f in source_files)
+    return newest
 def convert_type(source, dtype):
 …
     if dtype == F16:
         fbytes = 2
         source = _F16_PRAGMA + _convert_type(source, "half", "f")
+        source = _convert_type(source, "half", "f")
     elif dtype == F32:
         fbytes = 4
 …
     elif dtype == F64:
         fbytes = 8
         source = _F64_PRAGMA + source  # Source is already double
+        # no need to convert if it is already double
     elif dtype == F128:
         fbytes = 16
 …
+LOOP_OPEN = """\
+for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) {
+  const double %(name)s = loops[2*(%(name)s_i%(offset)s)];
+  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\
+_template_cache = {}
+def load_template(filename):
+    path = joinpath(TEMPLATE_ROOT, filename)
+    mtime = getmtime(path)
+    if filename not in _template_cache or mtime > _template_cache[filename][0]:
+        with open(path) as fid:
+            _template_cache[filename] = (mtime, fid.read(), path)
+    return _template_cache[filename][1]
+def model_templates():
+    # TODO: fails DRY; templates are listed in two places.
+    # should instead have model_info contain a list of paths
+    return [joinpath(TEMPLATE_ROOT, filename)
+            for filename in ('kernel_header.c', 'kernel_iq.c')]
+_FN_TEMPLATE = """\
+double %(name)s(%(pars)s);
+double %(name)s(%(pars)s) {
+    %(body)s
+}
 """
+def build_polydispersity_loops(pd_pars):
+    """
+    Build polydispersity loops
+    Returns loop opening and loop closing
+    """
+    depth = 4
+    offset = ""
+    loop_head = []
+    loop_end = []
+    for name in pd_pars:
+        subst = {'name': name, 'offset': offset}
+        loop_head.append(indent(LOOP_OPEN % subst, depth))
+        loop_end.insert(0, (" "*depth) + "}")
+        offset += '+N' + name
+        depth += 2
+    return "\n".join(loop_head), "\n".join(loop_end)
+C_KERNEL_TEMPLATE = None
+def _gen_fn(name, pars, body):
+    """
+    Generate a function given pars and body.
+    Returns the following string::
+         double fn(double a, double b, ...);
+         double fn(double a, double b, ...) {
+             ....
+         }
+    """
+    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
+    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl}
+def _call_pars(prefix, pars):
+    """
+    Return a list of *prefix.parameter* from parameter items.
+    """
+    return [p.as_call_reference(prefix) for p in pars]
+_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)",
+                           flags=re.MULTILINE)
+def _have_Iqxy(sources):
+    """
+    Return true if any file defines Iqxy.
+    Note this is not a C parser, and so can be easily confused by
+    non-standard syntax.  Also, it will incorrectly identify the following
+    as having Iqxy::
+        /*
+        double Iqxy(qx, qy, ...) { ... fill this in later ... }
+        */
+    If you want to comment out an Iqxy function, use // on the front of the
+    line instead.
+    """
+    for code in sources:
+        if _IQXY_PATTERN.search(code):
+            return True
+    else:
+        return False
 def make_source(model_info):
     """
 …
     # for computing volume even if we allow non-disperse volume parameters.
+    # Load template
+    global C_KERNEL_TEMPLATE
+    if C_KERNEL_TEMPLATE is None:
+        with open(C_KERNEL_TEMPLATE_PATH) as fid:
+            C_KERNEL_TEMPLATE = fid.read()
+    # Load additional sources
+    source = [open(f).read() for f in model_sources(model_info)]
+    # Prepare defines
+    defines = []
+    partype = model_info['partype']
+    pd_1d = partype['pd-1d']
+    pd_2d = partype['pd-2d']
+    fixed_1d = partype['fixed-1d']
+    fixed_2d = partype['fixed-1d']
+    iq_parameters = [p.name
+                     for p in model_info['parameters'][2:]  # skip scale, background
+                     if p.name in set(fixed_1d + pd_1d)]
+    iqxy_parameters = [p.name
+                       for p in model_info['parameters'][2:]  # skip scale, background
+                       if p.name in set(fixed_2d + pd_2d)]
+    volume_parameters = [p.name
+                         for p in model_info['parameters']
+                         if p.type == 'volume']
+    # Fill in defintions for volume parameters
+    if volume_parameters:
+        defines.append(('VOLUME_PARAMETERS',
+                        ','.join(volume_parameters)))
+        defines.append(('VOLUME_WEIGHT_PRODUCT',
+                        '*'.join(p + '_w' for p in volume_parameters)))
+    # Generate form_volume function from body only
+    partable = model_info['parameters']
+    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume.
+    # Note that scale and volume are not possible types.
+    # Load templates and user code
+    kernel_header = load_template('kernel_header.c')
+    kernel_code = load_template('kernel_iq.c')
+    user_code = [open(f).read() for f in model_sources(model_info)]
+    # Build initial sources
+    source = [kernel_header] + user_code
+    # Make parameters for q, qx, qy so that we can use them in declarations
+    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')]
+    # Generate form_volume function, etc. from body only
     if model_info['form_volume'] is not None:
+        if volume_parameters:
+            vol_par_decl = ', '.join('double ' + p for p in volume_parameters)
+        else:
+            vol_par_decl = 'void'
+        defines.append(('VOLUME_PARAMETER_DECLARATIONS',
+                        vol_par_decl))
+        fn = """\
+double form_volume(VOLUME_PARAMETER_DECLARATIONS);
+double form_volume(VOLUME_PARAMETER_DECLARATIONS) {
+    %(body)s
+}
+""" % {'body':model_info['form_volume']}
+        source.append(fn)
+    # Fill in definitions for Iq parameters
+    defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq'))
+    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters)))
+    if fixed_1d:
+        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS',
+                        ', \\\n    '.join('const double %s' % p for p in fixed_1d)))
+    if pd_1d:
+        defines.append(('IQ_WEIGHT_PRODUCT',
+                        '*'.join(p + '_w' for p in pd_1d)))
+        defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS',
+                        ', \\\n    '.join('const int N%s' % p for p in pd_1d)))
+        defines.append(('IQ_DISPERSION_LENGTH_SUM',
+                        '+'.join('N' + p for p in pd_1d)))
+        open_loops, close_loops = build_polydispersity_loops(pd_1d)
+        defines.append(('IQ_OPEN_LOOPS',
+                        open_loops.replace('\n', ' \\\n')))
+        defines.append(('IQ_CLOSE_LOOPS',
+                        close_loops.replace('\n', ' \\\n')))
+        pars = partable.form_volume_parameters
+        source.append(_gen_fn('form_volume', pars, model_info['form_volume']))
     if model_info['Iq'] is not None:
+        defines.append(('IQ_PARAMETER_DECLARATIONS',
+                        ', '.join('double ' + p for p in iq_parameters)))
+        fn = """\
+double Iq(double q, IQ_PARAMETER_DECLARATIONS);
+double Iq(double q, IQ_PARAMETER_DECLARATIONS) {
+    %(body)s
+}
+""" % {'body':model_info['Iq']}
+        source.append(fn)
+    # Fill in definitions for Iqxy parameters
+    defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy'))
+    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters)))
+    if fixed_2d:
+        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS',
+                        ', \\\n    '.join('const double %s' % p for p in fixed_2d)))
+    if pd_2d:
+        defines.append(('IQXY_WEIGHT_PRODUCT',
+                        '*'.join(p + '_w' for p in pd_2d)))
+        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS',
+                        ', \\\n    '.join('const int N%s' % p for p in pd_2d)))
+        defines.append(('IQXY_DISPERSION_LENGTH_SUM',
+                        '+'.join('N' + p for p in pd_2d)))
+        open_loops, close_loops = build_polydispersity_loops(pd_2d)
+        defines.append(('IQXY_OPEN_LOOPS',
+                        open_loops.replace('\n', ' \\\n')))
+        defines.append(('IQXY_CLOSE_LOOPS',
+                        close_loops.replace('\n', ' \\\n')))
+        pars = [q] + partable.iq_parameters
+        source.append(_gen_fn('Iq', pars, model_info['Iq']))
     if model_info['Iqxy'] is not None:
+        defines.append(('IQXY_PARAMETER_DECLARATIONS',
+                        ', '.join('double ' + p for p in iqxy_parameters)))
+        fn = """\
+double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS);
+double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) {
+    %(body)s
+}
+""" % {'body':model_info['Iqxy']}
+        source.append(fn)
+    # Need to know if we have a theta parameter for Iqxy; it is not there
+    # for the magnetic sphere model, for example, which has a magnetic
+    # orientation but no shape orientation.
+    if 'theta' in pd_2d:
+        defines.append(('IQXY_HAS_THETA', '1'))
+    #for d in defines: print(d)
+    defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines)
+    sources = '\n\n'.join(source)
+    return C_KERNEL_TEMPLATE % {
+        'DEFINES': defines,
+        'SOURCES': sources,
+        }
+def categorize_parameters(pars):
+    """
+    Build parameter categories out of the the parameter definitions.
+    Returns a dictionary of categories.
+    Note: these categories are subject to change, depending on the needs of
+    the UI and the needs of the kernel calling function.
+    The categories are as follows:
+    * *volume* list of volume parameter names
+    * *orientation* list of orientation parameters
+    * *magnetic* list of magnetic parameters
+    * *<empty string>* list of parameters that have no type info
+    Each parameter is in one and only one category.
+    The following derived categories are created:
+    * *fixed-1d* list of non-polydisperse parameters for 1D models
+    * *pd-1d* list of polydisperse parameters for 1D models
+    * *fixed-2d* list of non-polydisperse parameters for 2D models
+    * *pd-d2* list of polydisperse parameters for 2D models
+    """
+    partype = {
+        'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [],
+        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [],
+        'pd-rel': set(),
+    }
+    for p in pars:
+        if p.type == 'volume':
+            partype['pd-1d'].append(p.name)
+            partype['pd-2d'].append(p.name)
+            partype['pd-rel'].add(p.name)
+        elif p.type == 'magnetic':
+            partype['fixed-2d'].append(p.name)
+        elif p.type == 'orientation':
+            partype['pd-2d'].append(p.name)
+        elif p.type in ('', 'sld'):
+            partype['fixed-1d'].append(p.name)
+            partype['fixed-2d'].append(p.name)
+        else:
+            raise ValueError("unknown parameter type %r" % p.type)
+        partype[p.type].append(p.name)
+    return partype
+def process_parameters(model_info):
+    """
+    Process parameter block, precalculating parameter details.
+    """
+    # convert parameters into named tuples
+    for p in model_info['parameters']:
+        if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')):
+            p[4] = 'sld'
+            # TODO: make sure all models explicitly label their sld parameters
+            #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0]))
+    pars = [Parameter(*p) for p in model_info['parameters']]
+    # Fill in the derived attributes
+    model_info['parameters'] = pars
+    partype = categorize_parameters(pars)
+    model_info['limits'] = dict((p.name, p.limits) for p in pars)
+    model_info['partype'] = partype
+    model_info['defaults'] = dict((p.name, p.default) for p in pars)
+    if model_info.get('demo', None) is None:
+        model_info['demo'] = model_info['defaults']
+    model_info['has_2d'] = partype['orientation'] or partype['magnetic']
+        pars = [qx, qy] + partable.iqxy_parameters
+        source.append(_gen_fn('Iqxy', pars, model_info['Iqxy']))
+    # Define the parameter table
+    source.append("#define PARAMETER_TABLE \\")
+    source.append("\\\n".join(p.as_definition()
+                              for p in partable.kernel_parameters))
+    # Define the function calls
+    if partable.form_volume_parameters:
+        refs = _call_pars("_v.", partable.form_volume_parameters)
+        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs))
+    else:
+        # Model doesn't have volume.  We could make the kernel run a little
+        # faster by not using/transferring the volume normalizations, but
+        # the ifdef's reduce readability more than is worthwhile.
+        call_volume = "#define CALL_VOLUME(v) 1.0"
+    source.append(call_volume)
+    refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters)
+    call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs))
+    if _have_Iqxy(user_code):
+        # Call 2D model
+        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters)
+        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs))
+    else:
+        # Call 1D model with sqrt(qx^2 + qy^2)
+        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))")
+        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters)
+        pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:]
+        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt))
+    # Fill in definitions for numbers of parameters
+    source.append("#define MAX_PD %s"%partable.max_pd)
+    source.append("#define NPARS %d"%partable.npars)
+    # TODO: allow mixed python/opencl kernels?
+    # define the Iq kernel
+    source.append("#define KERNEL_NAME %s_Iq"%model_info['name'])
+    source.append(call_iq)
+    source.append(kernel_code)
+    source.append("#undef CALL_IQ")
+    source.append("#undef KERNEL_NAME")
+    # define the Iqxy kernel from the same source with different #defines
+    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name'])
+    source.append(call_iqxy)
+    source.append(kernel_code)
+    source.append("#undef CALL_IQ")
+    source.append("#undef KERNEL_NAME")
+    return '\n'.join(source)
+class CoordinationDetails(object):
+    def __init__(self, model_info):
+        parameters = model_info['parameters']
+        max_pd = parameters.max_pd
+        npars = parameters.npars
+        par_offset = 4*max_pd
+        self.details = np.zeros(par_offset + 3*npars + 4, 'i4')
+        # generate views on different parts of the array
+        self._pd_par     = self.details[0*max_pd:1*max_pd]
+        self._pd_length  = self.details[1*max_pd:2*max_pd]
+        self._pd_offset  = self.details[2*max_pd:3*max_pd]
+        self._pd_stride  = self.details[3*max_pd:4*max_pd]
+        self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars]
+        self._par_coord  = self.details[par_offset+1*npars:par_offset+2*npars]
+        self._pd_coord   = self.details[par_offset+2*npars:par_offset+3*npars]
+        # theta_par is fixed
+        self.details[-1] = parameters.theta_offset
+    @property
+    def ctypes(self): return self.details.ctypes
+    @property
+    def pd_par(self): return self._pd_par
+    @property
+    def pd_length(self): return self._pd_length
+    @property
+    def pd_offset(self): return self._pd_offset
+    @property
+    def pd_stride(self): return self._pd_stride
+    @property
+    def pd_coord(self): return self._pd_coord
+    @property
+    def par_coord(self): return self._par_coord
+    @property
+    def par_offset(self): return self._par_offset
+    @property
+    def num_active(self): return self.details[-4]
+    @num_active.setter
+    def num_active(self, v): self.details[-4] = v
+    @property
+    def total_pd(self): return self.details[-3]
+    @total_pd.setter
+    def total_pd(self, v): self.details[-3] = v
+    @property
+    def num_coord(self): return self.details[-2]
+    @num_coord.setter
+    def num_coord(self, v): self.details[-2] = v
+    @property
+    def theta_par(self): return self.details[-1]
+    def show(self):
+        print("total_pd", self.total_pd)
+        print("num_active", self.num_active)
+        print("pd_par", self.pd_par)
+        print("pd_length", self.pd_length)
+        print("pd_offset", self.pd_offset)
+        print("pd_stride", self.pd_stride)
+        print("par_offsets", self.par_offset)
+        print("num_coord", self.num_coord)
+        print("par_coord", self.par_coord)
+        print("pd_coord", self.pd_coord)
+        print("theta par", self.details[-1])
+def mono_details(model_info):
+    details = CoordinationDetails(model_info)
+    # The zero defaults for monodisperse systems are mostly fine
+    details.par_offset[:] = np.arange(2, len(details.par_offset)+2)
+    return details
+def poly_details(model_info, weights):
+    #print("weights",weights)
+    weights = weights[2:] # Skip scale and background
+    # Decreasing list of polydispersity lengths
+    # Note: the reversing view, x[::-1], does not require a copy
+    pd_length = np.array([len(w) for w in weights])
+    num_active = np.sum(pd_length>1)
+    if num_active > model_info['parameters'].max_pd:
+        raise ValueError("Too many polydisperse parameters")
+    pd_offset = np.cumsum(np.hstack((0, pd_length)))
+    idx = np.argsort(pd_length)[::-1][:num_active]
+    par_length = np.array([max(len(w),1) for w in weights])
+    pd_stride = np.cumprod(np.hstack((1, par_length[idx])))
+    par_offsets = np.cumsum(np.hstack((2, par_length)))
+    details = CoordinationDetails(model_info)
+    details.pd_par[:num_active] = idx
+    details.pd_length[:num_active] = pd_length[idx]
+    details.pd_offset[:num_active] = pd_offset[idx]
+    details.pd_stride[:num_active] = pd_stride[:-1]
+    details.par_offset[:] = par_offsets[:-1]
+    details.total_pd = pd_stride[-1]
+    details.num_active = num_active
+    # Without constraints coordinated parameters are just the pd parameters
+    details.par_coord[:num_active] = idx
+    details.pd_coord[:num_active] = 2**np.arange(num_active)
+    details.num_coord = num_active
+    #details.show()
+    return details
+def constrained_poly_details(model_info, weights, constraints):
+    # Need to find the independently varying pars and sort them
+    # Need to build a coordination list for the dependent variables
+    # Need to generate a constraints function which takes values
+    # and weights, returning par blocks
+    raise NotImplementedError("Can't handle constraints yet")
+def create_vector_Iq(model_info):
+    """
+    Define Iq as a vector function if it exists.
+    """
+    Iq = model_info['Iq']
+    if callable(Iq) and not getattr(Iq, 'vectorized', False):
+        def vector_Iq(q, *args):
+            """
+            Vectorized 1D kernel.
+            """
+            return np.array([Iq(qi, *args) for qi in q])
+        model_info['Iq'] = vector_Iq
+def create_vector_Iqxy(model_info):
+    """
+    Define Iqxy as a vector function if it exists, or default it from Iq().
+    """
+    Iq, Iqxy = model_info['Iq'], model_info['Iqxy']
+    if callable(Iqxy) and not getattr(Iqxy, 'vectorized', False):
+        def vector_Iqxy(qx, qy, *args):
+            """
+            Vectorized 2D kernel.
+            """
+            return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)])
+        model_info['Iqxy'] = vector_Iqxy
+    elif callable(Iq):
+        # Iq is vectorized because create_vector_Iq was already called.
+        def default_Iqxy(qx, qy, *args):
+            """
+            Default 2D kernel.
+            """
+            return Iq(np.sqrt(qx**2 + qy**2), *args)
+        model_info['Iqxy'] = default_Iqxy
+def create_default_functions(model_info):
+    """
+    Autogenerate missing functions, such as Iqxy from Iq.
+    This only works for Iqxy when Iq is written in python. :func:`make_source`
+    performs a similar role for Iq written in C.  This also vectorizes
+    any functions that are not already marked as vectorized.
+    """
+    create_vector_Iq(model_info)
+    create_vector_Iqxy(model_info)  # call create_vector_Iq() first
 def load_kernel_module(model_name):
 …
       model variants (e.g., the list of cases) or is None if there are no
       model variants
+    * *defaults* is the *{parameter: value}* table built from the parameter
+      description table.
+    * *limits* is the *{parameter: [min, max]}* table built from the
+      parameter description table.
+    * *partypes* categorizes the model parameters. See
+    * *par_type* categorizes the model parameters. See
       :func:`categorize_parameters` for details.
     * *demo* contains the *{parameter: value}* map used in compare (and maybe
 …
       *model_info* blocks for the composition objects.  This allows us to
       build complete product and mixture models from just the info.
     """
     # TODO: maybe turn model_info into a class ModelDefinition
+    parameters = COMMON_PARAMETERS + kernel_module.parameters
+    #print("make parameter table", kernel_module.parameters)
+    parameters = make_parameter_table(kernel_module.parameters)
     filename = abspath(kernel_module.__file__)
     kernel_id = splitext(basename(filename))[0]
 …
         filename=abspath(kernel_module.__file__),
         name=name,
         title=kernel_module.title,
         description=kernel_module.description,
+        title=getattr(kernel_module, 'title', name+" model"),
+        description=getattr(kernel_module, 'description', 'no description'),
         parameters=parameters,
         composition=None,
 …
         single=getattr(kernel_module, 'single', True),
         structure_factor=getattr(kernel_module, 'structure_factor', False),
+        profile_axes=getattr(kernel_module, 'profile_axes', ['x','y']),
         variant_info=getattr(kernel_module, 'invariant_info', None),
         demo=getattr(kernel_module, 'demo', None),
 …
         tests=getattr(kernel_module, 'tests', []),
+        )
     process_parameters(model_info)
+    set_demo(model_info, getattr(kernel_module, 'demo', None))
     # Check for optional functions
     functions = "ER VR form_volume Iq Iqxy shape sesans".split()
+    functions = "ER VR form_volume Iq Iqxy profile sesans".split()
     model_info.update((k, getattr(kernel_module, k, None)) for k in functions)
+    create_default_functions(model_info)
+    # Precalculate the monodisperse parameters
+    # TODO: make this a lazy evaluator
+    # make_model_info is called for every model on sasview startup
+    model_info['mono_details'] = mono_details(model_info)
     return model_info
 …
     Program which prints the source produced by the model.
     """
+    import sys
     if len(sys.argv) <= 1:
         print("usage: python -m sasmodels.generate modelname")

sasmodels/kernelcl.py

-                      r4d76711
+                      rd2bb604
 harmless, albeit annoying.
 """
+from __future__ import print_function
 import os
 import warnings
 …
 try:
+    raise NotImplementedError("OpenCL not yet implemented for new kernel template")
     import pyopencl as cl
     # Ask OpenCL for the default context so that we know that one exists
 …
 # of polydisperse parameters.
 MAX_LOOPS = 2048
+# Pragmas for enable OpenCL features.  Be sure to protect them so that they
+# still compile even if OpenCL is not present.
+_F16_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16)
+#  pragma OPENCL EXTENSION cl_khr_fp16: enable
+#endif
+"""
+_F64_PRAGMA = """\
+#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64)
+#  pragma OPENCL EXTENSION cl_khr_fp64: enable
+#endif
+"""
 …
         raise RuntimeError("%s not supported for devices"%dtype)
+    source = generate.convert_type(source, dtype)
+    source_list = [generate.convert_type(source, dtype)]
+    if dtype == generate.F16:
+        source_list.insert(0, _F16_PRAGMA)
+    elif dtype == generate.F64:
+        source_list.insert(0, _F64_PRAGMA)
     # Note: USE_SINCOS makes the intel cpu slower under opencl
     if context.devices[0].type == cl.device_type.GPU:
         source = "#define USE_SINCOS\n" + source
+        source_list.insert(0, "#define USE_SINCOS\n")
     options = (get_fast_inaccurate_build_options(context.devices[0])
                if fast else [])
+    source = "\n".join(source_list)
     program = cl.Program(context, source).build(options=options)
     return program
 …
         key = "%s-%s-%s"%(name, dtype, fast)
         if key not in self.compiled:
             #print("compiling",name)
+            print("compiling",name)
             dtype = np.dtype(dtype)
+            program = compile_model(self.get_context(dtype), source, dtype, fast)
+            program = compile_model(self.get_context(dtype),
+                                    str(source), dtype, fast)
             self.compiled[key] = program
         return self.compiled[key]
 …
         if self.program is None:
             compiler = environment().compile_program
             self.program = compiler(self.info['name'], self.source, self.dtype,
                                     self.fast)
+            self.program = compiler(self.info['name'], self.source,
+                                    self.dtype, self.fast)
         is_2d = len(q_vectors) == 2
         kernel_name = generate.kernel_name(self.info, is_2d)
 …
         # at this point, so instead using 32, which is good on the set of
         # architectures tested so far.
+        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
+        if self.is_2d:
+            # Note: 17 rather than 15 because results is 2 elements
+            # longer than input.
+            width = ((self.nq+17)//16)*16
+            self.q = np.empty((width, 2), dtype=dtype)
+            self.q[:self.nq, 0] = q_vectors[0]
+            self.q[:self.nq, 1] = q_vectors[1]
+        else:
+            # Note: 33 rather than 31 because results is 2 elements
+            # longer than input.
+            width = ((self.nq+33)//32)*32
+            self.q = np.empty(width, dtype=dtype)
+            self.q[:self.nq] = q_vectors[0]
+        self.global_size = [self.q.shape[0]]
         context = env.get_context(self.dtype)
-        self.global_size = [self.q_vectors[0].size]
         #print("creating inputs of size", self.global_size)
+        self.q_buffers = [
+            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
+            for q in self.q_vectors
+        ]
+        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                             hostbuf=self.q)
     def release(self):
 …
         Free the memory.
         """
         for b in self.q_buffers:
             b.release()
         self.q_buffers = []
+        if self.q is not None:
+            self.q.release()
+            self.q = None
     def __del__(self):
 …
     """
     def __init__(self, kernel, model_info, q_vectors, dtype):
+        max_pd = model_info['max_pd']
+        npars = len(model_info['parameters'])-2
         q_input = GpuInput(q_vectors, dtype)
+        self.dtype = dtype
+        self.dim = '2d' if q_input.is_2d else '1d'
         self.kernel = kernel
         self.info = model_info
+        self.res = np.empty(q_input.nq, q_input.dtype)
+        dim = '2d' if q_input.is_2d else '1d'
+        self.fixed_pars = model_info['partype']['fixed-' + dim]
+        self.pd_pars = model_info['partype']['pd-' + dim]
+        self.pd_stop_index = 4*max_pd-1
+        # plus three for the normalization values
+        self.result = np.empty(q_input.nq+3, q_input.dtype)
         # Inputs and outputs for each kernel call
 …
         env = environment()
         self.queue = env.get_queue(dtype)
+        self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
+* MAX_LOOPS * q_input.dtype.itemsize)
+        self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
+        # details is int32 data, padded to an 8 integer boundary
+        size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8
+        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
                                q_input.global_size[0] * q_input.dtype.itemsize)
         self.q_input = q_input
         self._need_release = [self.loops_b, self.res_b, self.q_input]
     def __call__(self, fixed_pars, pd_pars, cutoff):
+        self.q_input = q_input # allocated by GpuInput above
+        self._need_release = [ self.result_b, self.q_input ]
+    def __call__(self, details, weights, values, cutoff):
         real = (np.float32 if self.q_input.dtype == generate.F32
                 else np.float64 if self.q_input.dtype == generate.F64
                 else np.float16 if self.q_input.dtype == generate.F16
                 else np.float32)  # will never get here, so use np.float32
+        #print "pars", fixed_pars, pd_pars
+        res_bi = self.res_b
+        nq = np.uint32(self.q_input.nq)
+        if pd_pars:
+            cutoff = real(cutoff)
+            loops_N = [np.uint32(len(p[0])) for p in pd_pars]
+            loops = np.hstack(pd_pars) \
+                if pd_pars else np.empty(0, dtype=self.q_input.dtype)
+            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten()
+            #print("loops",Nloops, loops)
+            #import sys; print("opencl eval",pars)
+            #print("opencl eval",pars)
+            if len(loops) > 2 * MAX_LOOPS:
+                raise ValueError("too many polydispersity points")
+            loops_bi = self.loops_b
+            cl.enqueue_copy(self.queue, loops_bi, loops)
+            loops_l = cl.LocalMemory(len(loops.data))
+            #ctx = environment().context
+            #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops)
+            dispersed = [loops_bi, loops_l, cutoff] + loops_N
+        else:
+            dispersed = []
+        fixed = [real(p) for p in fixed_pars]
+        args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed
+        assert details.dtype == np.int32
+        assert weights.dtype == real and values.dtype == real
+        context = self.queue.context
+        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                              hostbuf=details)
+        weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                              hostbuf=weights)
+        values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                             hostbuf=values)
+        start, stop = 0, self.details[self.pd_stop_index]
+        args = [
+            np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop),
+            self.details_b, self.weights_b, self.values_b,
+            self.q_input.q_b, self.result_b, real(cutoff),
+        ]
         self.kernel(self.queue, self.q_input.global_size, None, *args)
+        cl.enqueue_copy(self.queue, self.res, res_bi)
+        return self.res
+        cl.enqueue_copy(self.queue, self.result, self.result_b)
+        [v.release() for v in details_b, weights_b, values_b]
+        return self.result[:self.nq]
     def release(self):

sasmodels/kerneldll.py

-                      r4d76711
+                      r1e2a1ba
 import tempfile
 import ctypes as ct
+from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float
+import _ctypes
+from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float
 import numpy as np
 …
         COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
         if "SAS_OPENMP" in os.environ:
             COMPILE = COMPILE + " -fopenmp"
+            COMPILE += " -fopenmp"
 else:
     COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
 …
 def dll_path(model_info, dtype="double"):
     """
     Path to the compiled model defined by *model_info*.
     """
     from os.path import join as joinpath, split as splitpath, splitext
     basename = splitext(splitpath(model_info['filename'])[1])[0]
     if np.dtype(dtype) == generate.F32:
+        basename += "32"
     elif np.dtype(dtype) == generate.F64:
         basename += "64"
     else:
         basename += "128"
     return joinpath(DLL_PATH, basename+'.so')
+def dll_name(model_info, dtype):
+    """
+    Name of the dll containing the model.  This is the base file name without
+    any path or extension, with a form such as 'sas_sphere32'.
+    """
+    bits = 8*dtype.itemsize
+    return "sas_%s%d"%(model_info['id'], bits)
+def dll_path(model_info, dtype):
+    """
+    Complete path to the dll for the model.  Note that the dll may not
+    exist yet if it hasn't been compiled.
+    """
+    return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so")
 def make_dll(source, model_info, dtype="double"):
 …
         dtype = generate.F64  # Force 64-bit dll
-    if dtype == generate.F32: # 32-bit dll
-        tempfile_prefix = 'sas_' + model_info['name'] + '32_'
-    elif dtype == generate.F64:
-        tempfile_prefix = 'sas_' + model_info['name'] + '64_'
-    else:
-        tempfile_prefix = 'sas_' + model_info['name'] + '128_'
     source = generate.convert_type(source, dtype)
     source_files = generate.model_sources(model_info) + [model_info['filename']]
+    newest = generate.timestamp(model_info)
     dll = dll_path(model_info, dtype)
-    newest = max(os.path.getmtime(f) for f in source_files)
     if not os.path.exists(dll) or os.path.getmtime(dll) < newest:
         # Replace with a proper temp file
         fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix)
+        basename = dll_name(model_info, dtype) + "_"
+        fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename)
         os.fdopen(fid, "w").write(source)
         command = COMPILE%{"source":filename, "output":dll}
 …
     return DllModel(filename, model_info, dtype=dtype)
-IQ_ARGS = [c_void_p, c_void_p, c_int]
-IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int]
 class DllModel(object):
     """
 …
     def _load_dll(self):
-        Nfixed1d = len(self.info['partype']['fixed-1d'])
-        Nfixed2d = len(self.info['partype']['fixed-2d'])
-        Npd1d = len(self.info['partype']['pd-1d'])
-        Npd2d = len(self.info['partype']['pd-2d'])
         #print("dll", self.dllpath)
         try:
 …
               else c_double if self.dtype == generate.F64
               else c_longdouble)
+        pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else []
+        pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []
+        # int, int, int, int*, double*, double*, double*, double*, double*, double
+        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp]
         self.Iq = self.dll[generate.kernel_name(self.info, False)]
-        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d
         self.Iqxy = self.dll[generate.kernel_name(self.info, True)]
+        self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d
+        self.release()
+        self.Iq.argtypes = argtypes
+        self.Iqxy.argtypes = argtypes
     def __getstate__(self):
 …
     """
     def __init__(self, kernel, model_info, q_input):
+        self.kernel = kernel
         self.info = model_info
         self.q_input = q_input
+        self.kernel = kernel
+        self.res = np.empty(q_input.nq, q_input.dtype)
+        dim = '2d' if q_input.is_2d else '1d'
+        self.fixed_pars = model_info['partype']['fixed-' + dim]
+        self.pd_pars = model_info['partype']['pd-' + dim]
+        # In dll kernel, but not in opencl kernel
+        self.p_res = self.res.ctypes.data
+    def __call__(self, fixed_pars, pd_pars, cutoff):
+        self.dtype = q_input.dtype
+        self.dim = '2d' if q_input.is_2d else '1d'
+        self.result = np.empty(q_input.nq+3, q_input.dtype)
+    def __call__(self, details, weights, values, cutoff):
         real = (np.float32 if self.q_input.dtype == generate.F32
                 else np.float64 if self.q_input.dtype == generate.F64
                 else np.float128)
+        nq = c_int(self.q_input.nq)
+        if pd_pars:
+            cutoff = real(cutoff)
+            loops_N = [np.uint32(len(p[0])) for p in pd_pars]
+            loops = np.hstack(pd_pars)
+            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten()
+            p_loops = loops.ctypes.data
+            dispersed = [p_loops, cutoff] + loops_N
+        else:
+            dispersed = []
+        fixed = [real(p) for p in fixed_pars]
+        args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed
+        #print(pars)
+        assert isinstance(details, generate.CoordinationDetails)
+        assert weights.dtype == real and values.dtype == real
+        start, stop = 0, details.total_pd
+        #print("in kerneldll")
+        #print("weights", weights)
+        #print("values", values)
+        args = [
+            self.q_input.nq, # nq
+            start, # pd_start
+            stop, # pd_stop pd_stride[MAX_PD]
+            details.ctypes.data, # problem
+            weights.ctypes.data,  # weights
+            values.ctypes.data,  #pars
+            self.q_input.q.ctypes.data, #q
+            self.result.ctypes.data,   # results
+            real(cutoff), # cutoff
+            ]
         self.kernel(*args)
+        return self.res
+        return self.result[:-3]
     def release(self):

sasmodels/kernelpy.py

-                      r4d76711
+                      r6e7ff6d
         self.dtype = dtype
         self.is_2d = (len(q_vectors) == 2)
+        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors]
+        self.q_pointers = [q.ctypes.data for q in self.q_vectors]
+        if self.is_2d:
+            self.q = np.empty((self.nq, 2), dtype=dtype)
+            self.q[:, 0] = q_vectors[0]
+            self.q[:, 1] = q_vectors[1]
+        else:
+            self.q = np.empty(self.nq, dtype=dtype)
+            self.q[:self.nq] = q_vectors[0]
     def release(self):
 …
         Free resources associated with the model inputs.
         """
         self.q_vectors = []
+        self.q = None
 class PyKernel(object):
 …
     """
     def __init__(self, kernel, model_info, q_input):
+        self.dtype = np.dtype('d')
         self.info = model_info
         self.q_input = q_input
         self.res = np.empty(q_input.nq, q_input.dtype)
+        dim = '2d' if q_input.is_2d else '1d'
+        # Loop over q unless user promises that the kernel is vectorized by
+        # taggining it with vectorized=True
+        if not getattr(kernel, 'vectorized', False):
+            if dim == '2d':
+                def vector_kernel(qx, qy, *args):
+                    """
+                    Vectorized 2D kernel.
+                    """
+                    return np.array([kernel(qxi, qyi, *args)
+                                     for qxi, qyi in zip(qx, qy)])
+        self.kernel = kernel
+        self.dim = '2d' if q_input.is_2d else '1d'
+        partable = model_info['parameters']
+        kernel_parameters = (partable.iqxy_parameters if q_input.is_2d
+                             else partable.iq_parameters)
+        volume_parameters = partable.form_volume_parameters
+        # Create an array to hold the parameter values.  There will be a
+        # single array whose values are updated as the calculator goes
+        # through the loop.  Arguments to the kernel and volume functions
+        # will use views into this vector, relying on the fact that a
+        # an array of no dimensions acts like a scalar.
+        parameter_vector = np.empty(len(partable.call_parameters), 'd')
+        # Create views into the array to hold the arguments
+        offset = 2
+        kernel_args, volume_args = [], []
+        for p in partable.kernel_parameters:
+            if p.length == 1:
+                # Scalar values are length 1 vectors with no dimensions.
+                v = parameter_vector[offset:offset+1].reshape(())
             else:
+                def vector_kernel(q, *args):
+                    """
+                    Vectorized 1D kernel.
+                    """
+                    return np.array([kernel(qi, *args)
+                                     for qi in q])
+            self.kernel = vector_kernel
+                # Vector values are simple views.
+                v = parameter_vector[offset:offset+p.length]
+            offset += p.length
+            if p in kernel_parameters:
+                kernel_args.append(v)
+            if p in volume_parameters:
+                volume_args.append(p)
+        # Hold on to the parameter vector so we can use it to call kernel later.
+        # This may also be required to preserve the views into the vector.
+        self._parameter_vector = parameter_vector
+        # Generate a closure which calls the kernel with the views into the
+        # parameter array.
+        if q_input.is_2d:
+            form = model_info['Iqxy']
+            qx, qy = q_input.q[:,0], q_input.q[:,1]
+            self._form = lambda: form(qx, qy, *kernel_args)
         else:
+            self.kernel = kernel
+        fixed_pars = model_info['partype']['fixed-' + dim]
+        pd_pars = model_info['partype']['pd-' + dim]
+        vol_pars = model_info['partype']['volume']
+        # First two fixed pars are scale and background
+        pars = [p.name for p in model_info['parameters'][2:]]
+        offset = len(self.q_input.q_vectors)
+        self.args = self.q_input.q_vectors + [None] * len(pars)
+        self.fixed_index = np.array([pars.index(p) + offset
+                                     for p in fixed_pars[2:]])
+        self.pd_index = np.array([pars.index(p) + offset
+                                  for p in pd_pars])
+        self.vol_index = np.array([pars.index(p) + offset
+                                   for p in vol_pars])
+        try: self.theta_index = pars.index('theta') + offset
+        except ValueError: self.theta_index = -1
+        # Caller needs fixed_pars and pd_pars
+        self.fixed_pars = fixed_pars
+        self.pd_pars = pd_pars
+    def __call__(self, fixed, pd, cutoff=1e-5):
+        #print("fixed",fixed)
+        #print("pd", pd)
+        args = self.args[:]  # grab a copy of the args
+        form, form_volume = self.kernel, self.info['form_volume']
+        # First two fixed
+        scale, background = fixed[:2]
+        for index, value in zip(self.fixed_index, fixed[2:]):
+            args[index] = float(value)
+        res = _loops(form, form_volume, cutoff, scale, background, args,
+                     pd, self.pd_index, self.vol_index, self.theta_index)
+            form = model_info['Iq']
+            q = q_input.q
+            self._form = lambda: form(q, *kernel_args)
+        # Generate a closure which calls the form_volume if it exists.
+        form_volume = model_info['form_volume']
+        self._volume = ((lambda: form_volume(*volume_args)) if form_volume
+                        else (lambda: 1.0))
+    def __call__(self, details, weights, values, cutoff):
+        # type: (.generate.CoordinationDetails, np.ndarray, np.ndarray, float) -> np.ndarray
+        res = _loops(self._parameter_vector, self._form, self._volume,
+                     self.q_input.nq, details, weights, values, cutoff)
         return res
 …
         self.q_input = None
+def _loops(form, form_volume, cutoff, scale, background,
+           args, pd, pd_index, vol_index, theta_index):
+    """
+    *form* is the name of the form function, which should be vectorized over
+    q, but otherwise have an interface like the opencl kernels, with the
+    q parameters first followed by the individual parameters in the order
+    given in model.parameters (see :mod:`sasmodels.generate`).
+    *form_volume* calculates the volume of the shape.  *vol_index* gives
+    the list of volume parameters
+    *cutoff* ignores the corners of the dispersion hypercube
+    *scale*, *background* multiplies the resulting form and adds an offset
+    *args* is the prepopulated set of arguments to the form function, starting
+    with the q vectors, and including slots for all the parameters.  The
+    values for the parameters will be substituted with values from the
+    dispersion functions.
+    *pd* is the list of dispersion parameters
+    *pd_index* are the indices of the dispersion parameters in *args*
+    *vol_index* are the indices of the volume parameters in *args*
+    *theta_index* is the index of the theta parameter for the sasview
+    spherical correction, or -1 if there is no angular dispersion
+    """
+def _loops(parameters,  # type: np.ndarray
+           form,        # type: Callable[[], np.ndarray]
+           form_volume, # type: Callable[[], float]
+           nq,          # type: int
+           details,     # type: .generate.CoordinationDetails
+           weights,     # type: np.ndarray
+           values,      # type: np.ndarray
+           cutoff,      # type: float
+           ):           # type: (...) -> None
     ################################################################
     #                                                              #
 …
     #                                                              #
     ################################################################
+    weight = np.empty(len(pd), 'd')
+    if weight.size > 0:
+        # weight vector, to be populated by polydispersity loops
+        # identify which pd parameters are volume parameters
+        vol_weight_index = np.array([(index in vol_index) for index in pd_index])
+        # Sort parameters in decreasing order of pd length
+        Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
+        order = np.argsort(Npd)[::-1]
+        stride = np.cumprod(Npd[order])
+        pd = [pd[index] for index in order]
+        pd_index = pd_index[order]
+        vol_weight_index = vol_weight_index[order]
+        fast_value = pd[0][0]
+        fast_weight = pd[0][1]
+    parameters[2:] = values[details.par_offset]
+    scale, background = values[0], values[1]
+    if details.num_active == 0:
+        norm = float(form_volume())
+        if norm > 0.0:
+            return (scale/norm)*form() + background
+        else:
+            return np.ones(nq, 'd')*background
+    partial_weight = np.NaN
+    spherical_correction = 1.0
+    pd_stride = details.pd_stride[:details.num_active]
+    pd_length = details.pd_length[:details.num_active]
+    pd_offset = details.pd_offset[:details.num_active]
+    pd_index = np.empty_like(pd_offset)
+    offset = np.empty_like(details.par_offset)
+    theta = details.theta_par
+    fast_length = pd_length[0]
+    pd_index[0] = fast_length
+    total = np.zeros(nq, 'd')
+    norm = 0.0
+    for loop_index in range(details.total_pd):
+        # update polydispersity parameter values
+        if pd_index[0] == fast_length:
+            pd_index[:] = (loop_index/pd_stride)%pd_length
+            partial_weight = np.prod(weights[pd_offset+pd_index][1:])
+            for k in range(details.num_coord):
+                par = details.par_coord[k]
+                coord = details.pd_coord[k]
+                this_offset = details.par_offset[par]
+                block_size = 1
+                for bit in xrange(32):
+                    if coord&1:
+                        this_offset += block_size * pd_index[bit]
+                        block_size *= pd_length[bit]
+                    coord >>= 1
+                    if coord == 0: break
+                offset[par] = this_offset
+                parameters[par] = values[this_offset]
+                if par == theta and not (details.par_coord[k]&1):
+                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6)
+        for k in range(details.num_coord):
+            if details.pd_coord[k]&1:
+                par = details.par_coord[k]
+                parameters[par] = values[offset[par]]
+                offset[par] += 1
+                if par == theta:
+                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6)
+        weight = partial_weight * weights[pd_offset[0] + pd_index[0]]
+        pd_index[0] += 1
+        if weight > cutoff:
+            # Call the scattering function
+            # Assume that NaNs are only generated if the parameters are bad;
+            # exclude all q for that NaN.  Even better would be to have an
+            # INVALID expression like the C models, but that is too expensive.
+            I = form()
+            if np.isnan(I).any(): continue
+            # update value and norm
+            weight *= spherical_correction
+            total += weight * I
+            norm += weight * form_volume()
+    if norm > 0.0:
+        return (scale/norm)*total + background
     else:
+        stride = np.array([1])
+        vol_weight_index = slice(None, None)
+        # keep lint happy
+        fast_value = [None]
+        fast_weight = [None]
+    ret = np.zeros_like(args[0])
+    norm = np.zeros_like(ret)
+    vol = np.zeros_like(ret)
+    vol_norm = np.zeros_like(ret)
+    for k in range(stride[-1]):
+        # update polydispersity parameter values
+        fast_index = k % stride[0]
+        if fast_index == 0:  # bottom loop complete ... check all other loops
+            if weight.size > 0:
+                for i, index, in enumerate(k % stride):
+                    args[pd_index[i]] = pd[i][0][index]
+                    weight[i] = pd[i][1][index]
+        else:
+            args[pd_index[0]] = fast_value[fast_index]
+            weight[0] = fast_weight[fast_index]
+        # Computes the weight, and if it is not sufficient then ignore this
+        # parameter set.
+        # Note: could precompute w1*...*wn so we only need to multiply by w0
+        w = np.prod(weight)
+        if w > cutoff:
+            # Note: can precompute spherical correction if theta_index is not
+            # the fast index. Correction factor for spherical integration
+            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0
+            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2
+                                    if theta_index >= 0 else 1.0)
+            #spherical_correction = 1.0
+            # Call the scattering function and adds it to the total.
+            I = form(*args)
+            if np.isnan(I).any(): continue
+            ret += w * I * spherical_correction
+            norm += w
+            # Volume normalization.
+            # If there are "volume" polydispersity parameters, then these
+            # will be used to call the form_volume function from the user
+            # supplied kernel, and accumulate a normalized weight.
+            # Note: can precompute volume norm if fast index is not a volume
+            if form_volume:
+                vol_args = [args[index] for index in vol_index]
+                vol_weight = np.prod(weight[vol_weight_index])
+                vol += vol_weight * form_volume(*vol_args)
+                vol_norm += vol_weight
+    positive = (vol * vol_norm != 0.0)
+    ret[positive] *= vol_norm[positive] / vol[positive]
+    result = scale * ret / norm + background
+    return result
+        return np.ones(nq, 'd')*background

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
 import numpy as np
 from .generate import process_parameters, COMMON_PARAMETERS, Parameter
+from .modelinfo import Parameter, ParameterTable
 SCALE=0
 …
     # Build new parameter list
     pars = COMMON_PARAMETERS + []
+    pars = []
     for k, part in enumerate(parts):
         # Parameter prefix per model, A_, B_, ...
+        # Note that prefix must also be applied to id and length_control
+        # to support vector parameters
         prefix = chr(ord('A')+k) + '_'
+        for p in part['parameters']:
+            # No background on the individual mixture elements
+            if p.name == 'background':
+                continue
+            # TODO: promote Parameter to a full class
+            # this code knows too much about the implementation!
+            p = list(p)
+            if p[0] == 'scale':  # allow model subtraction
+                p[3] = [-np.inf, np.inf]  # limits
+            p[0] = prefix+p[0]   # relabel parameters with A_, ...
+        pars.append(Parameter(prefix+'scale'))
+        for p in part['parameters'].kernel_pars:
+            p = copy(p)
+            p.name = prefix+p.name
+            p.id = prefix+p.id
+            if p.length_control is not None:
+                p.length_control = prefix+p.length_control
             pars.append(p)
+    partable = ParameterTable(pars)
     model_info = {}
 …
     model_info['docs'] = model_info['title']
     model_info['category'] = "custom"
     model_info['parameters'] = pars
+    model_info['parameters'] = partable
     #model_info['single'] = any(part['single'] for part in parts)
     model_info['structure_factor'] = False
 …
     # Remember the component info blocks so we can build the model
     model_info['composition'] = ('mixture', parts)
-    process_parameters(model_info)
-    return model_info

sasmodels/models/cylinder.c

-                      r26141cb
+                      re9b1663d
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
     double radius, double length, double theta, double phi);
+#define INVALID(v) (v.radius<0 || v.length<0)
 double form_volume(double radius, double length)
 …
     double length)
+{
-    // TODO: return NaN if radius<0 or length<0?
     // precompute qr and qh to save time in the loop
     const double qr = q*radius;
 …
     double phi)
+{
-    // TODO: return NaN if radius<0 or length<0?
     double sn, cn; // slots to hold sincos function output

sasmodels/models/flexible_cylinder.c

-                      re6408d0
+                      r4937980
+double form_volume(double length, double kuhn_length, double radius);
+double Iq(double q, double length, double kuhn_length, double radius,
+          double sld, double solvent_sld);
+double Iqxy(double qx, double qy, double length, double kuhn_length,
+            double radius, double sld, double solvent_sld);
+double flexible_cylinder_kernel(double q, double length, double kuhn_length,
+                                double radius, double sld, double solvent_sld);
+double form_volume(double length, double kuhn_length, double radius)
+static double
+form_volume(length, kuhn_length, radius)
+{
     return 1.0;
+}
+double flexible_cylinder_kernel(double q,
+          double length,
+          double kuhn_length,
+          double radius,
+          double sld,
+          double solvent_sld)
+static double
+Iq(double q,
+   double length,
+   double kuhn_length,
+   double radius,
+   double sld,
+   double solvent_sld)
+{
+    const double cont = sld-solvent_sld;
+    const double qr = q*radius;
+    //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr);
+    const double crossSect = sas_J1c(qr);
+    double flex = Sk_WR(q,length,kuhn_length);
+    flex *= crossSect*crossSect;
+    flex *= M_PI*radius*radius*length;
+    flex *= cont*cont;
+    flex *= 1.0e-4;
+    return flex;
+    const double contrast = sld - solvent_sld;
+    const double cross_section = sas_J1c(q*radius);
+    const double volume = M_PI*radius*radius*length;
+    const double flex = Sk_WR(q, length, kuhn_length);
+    return 1.0e-4 * volume * square(contrast*cross_section) * flex;
+}
-double Iq(double q,
-          double length,
-          double kuhn_length,
-          double radius,
-          double sld,
-          double solvent_sld)
+{
-    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
-    return result;
+}
-double Iqxy(double qx, double qy,
-            double length,
-            double kuhn_length,
-            double radius,
-            double sld,
-            double solvent_sld)
+{
-    double q;
-    q = sqrt(qx*qx+qy*qy);
-    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
-    return result;
+}

sasmodels/models/gel_fit.c

-                      r30b4ddf
+                      r03cac08
+double form_volume(void);
+double Iq(double q,
+          double guinier_scale,
+          double lorentzian_scale,
+          double gyration_radius,
+          double fractal_exp,
+          double cor_length);
+double Iqxy(double qx, double qy,
+          double guinier_scale,
+          double lorentzian_scale,
+          double gyration_radius,
+          double fractal_exp,
+          double cor_length);
+static double _gel_fit_kernel(double q,
+static double Iq(double q,
           double guinier_scale,
           double lorentzian_scale,
 …
     // Lorentzian Term
     ////////////////////////double a(x[i]*x[i]*zeta*zeta);
     double lorentzian_term = q*q*cor_length*cor_length;
+    double lorentzian_term = square(q*cor_length);
     lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term;
     lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) );
 …
     // Exponential Term
     ////////////////////////double d(x[i]*x[i]*rg*rg);
     double exp_term = q*q*gyration_radius*gyration_radius;
+    double exp_term = square(q*gyration_radius);
     exp_term = exp(-1.0*(exp_term/3.0));
 …
     return result;
+}
-double form_volume(void){
-    // Unused, so free to return garbage.
-    return NAN;
+}
-double Iq(double q,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
+{
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
+}
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-          double guinier_scale,
-          double lorentzian_scale,
-          double gyration_radius,
-          double fractal_exp,
-          double cor_length)
+{
-    double q = sqrt(qx*qx + qy*qy);
-    return _gel_fit_kernel(q,
-                          guinier_scale,
-                          lorentzian_scale,
-                          gyration_radius,
-                          fractal_exp,
-                          cor_length);
+}

sasmodels/models/hardsphere.py

-                      rec45c4f
+                      rd2bb604
    """
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/hayter_msa.py

-                      rec45c4f
+                      rd2bb604
     return 1.0;
     """
-Iqxy = """
-    // never called since no orientation or magnetic parameters.
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lamellar.py

-                      rec45c4f
+                      rd2bb604
     """
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lamellar_hg.py

-                      rec45c4f
+                      rd2bb604
     """
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lamellar_hg_stack_caille.py

-                      rec45c4f
+                      rd2bb604
     """
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lamellar_stack_caille.py

-                      rec45c4f
+                      rd2bb604
     """
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lamellar_stack_paracrystal.py

-                      rec45c4f
+                      rd2bb604
     """
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lib/sas_JN.c

-                      re6408d0
+                      r4937980
 */
+static double
+sas_JN( int n, double x ) {
+double sas_JN( int n, double x );
+double sas_JN( int n, double x ) {
     const double MACHEP = 1.11022302462515654042E-16;

sasmodels/models/lib/sph_j1c.c

-                      re6f1410
+                      rba32cdd
 * using double precision that are the source.
 */
+double sph_j1c(double q);
 // The choice of the number of terms in the series and the cutoff value for
 …
 #endif
-double sph_j1c(double q);
 double sph_j1c(double q)
+{

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
     Functions for WRC implementation of flexible cylinders
 */
+double Sk_WR(double q, double L, double b);
 static double
 AlphaSquare(double x)
 …
+}
-double Sk_WR(double q, double L, double b);
 double Sk_WR(double q, double L, double b)
+{

sasmodels/models/onion.c

-                      rfdb1487
+                      rce896fd
     double thickness, double A)
+{
   const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double alpha = A * r/thickness;
 …
     double sinqr, cosqr;
     SINCOS(qr, sinqr, cosqr);
     fun = -3.0(
+    fun = -3.0*(
             ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
                 - (alpha*sinqr/qr - cosqr)
 …
 f_linear(double q, double r, double sld, double slope)
+{
   const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double bes = sph_j1c(qr);
 …
+{
   const double bes = sph_j1c(q * r);
   const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   return sld * vol * bes;
+}
 …
     r += thickness[i];
+  }
   return 4.0/3.0 * M_PI * r * r * r;
+  return M_4PI_3*cube(r);
+}
 static double
 Iq(double q, double core_sld, double core_radius, double solvent_sld,
     double n, double in_sld[], double out_sld[], double thickness[],
+Iq(double q, double sld_core, double core_radius, double sld_solvent,
+    double n, double sld_in[], double sld_out[], double thickness[],
     double A[])
+{
   int i;
   r = core_radius;
   f = f_constant(q, r, core_sld);
+  double r = core_radius;
+  double f = f_constant(q, r, sld_core);
   for (i=0; i<n; i++){
     const double r0 = r;
 …
+    }
+  }
   f -= f_constant(q, r, solvent_sld);
   f2 = f * f * 1.0e-4;
+  f -= f_constant(q, r, sld_solvent);
+  const double f2 = f * f * 1.0e-4;
   return f2;

sasmodels/models/onion.py

-                      rec45c4f
+                      rea05c87
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
                "Core scattering length density"],
               ["core_radius", "Ang", 200., [0, inf], "volume",
                "Radius of the core"],
               ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "",
+              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
                "Solvent scattering length density"],
               ["n", "", 1, [0, 10], "volume",
                "number of shells"],
               ["in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+              ["sld_in[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
                "scattering length density at the inner radius of shell k"],
               ["out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
+              ["sld_out[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "",
                "scattering length density at the outer radius of shell k"],
               ["thickness[n]", "Ang", 40., [0, inf], "volume",
 …
+              ]
+#source = ["lib/sph_j1c.c", "onion.c"]
+def Iq(q, *args, **kw):
+    return q
+def Iqxy(qx, *args, **kw):
+    return qx
+def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A):
+source = ["lib/sph_j1c.c", "onion.c"]
+#def Iq(q, *args, **kw):
+#    return q
+profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
+def profile(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A):
     """
     SLD profile
 …
 demo = {
     "solvent_sld": 2.2,
     "core_sld": 1.0,
+    "sld_solvent": 2.2,
+    "sld_core": 1.0,
     "core_radius": 100,
     "n": 4,
     "in_sld": [0.5, 1.5, 0.9, 2.0],
     "out_sld": [nan, 0.9, 1.2, 1.6],
+    "sld_in": [0.5, 1.5, 0.9, 2.0],
+    "sld_out": [nan, 0.9, 1.2, 1.6],
     "thickness": [50, 75, 150, 75],
     "A": [0, -1, 1e-4, 1],

sasmodels/models/rpa.c

-                      r13ed84c
+                      rd2bb604
 double Iq(double q, double case_num,
+    double Na, double Phia, double va, double a_sld, double ba,
+    double Nb, double Phib, double vb, double b_sld, double bb,
+    double Nc, double Phic, double vc, double c_sld, double bc,
+    double Nd, double Phid, double vd, double d_sld, double bd,
+    double N[], double Phi[], double v[], double L[], double b[],
     double Kab, double Kac, double Kad,
     double Kbc, double Kbd, double Kcd
     );
-double Iqxy(double qx, double qy, double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
-    double Kab, double Kac, double Kad,
-    double Kbc, double Kbd, double Kcd
-    );
-double form_volume(void);
-double form_volume(void)
+{
-    return 1.0;
+}
 double Iq(double q, double case_num,
+    double Na, double Phia, double va, double La, double ba,
+    double Nb, double Phib, double vb, double Lb, double bb,
+    double Nc, double Phic, double vc, double Lc, double bc,
+    double Nd, double Phid, double vd, double Ld, double bd,
+    double N[], double Phi[], double v[], double L[], double b[],
     double Kab, double Kac, double Kad,
     double Kbc, double Kbd, double Kcd
 …
 #if 0  // Sasview defaults
   if (icase <= 1) {
     Na=Nb=1000.0;
     Phia=Phib=0.0000001;
+    N[0]=N[1]=1000.0;
+    Phi[0]=Phi[1]=0.0000001;
     Kab=Kac=Kad=Kbc=Kbd=-0.0004;
     La=Lb=1e-12;
     va=vb=100.0;
     ba=bb=5.0;
+    L[0]=L[1]=1e-12;
+    v[0]=v[1]=100.0;
+    b[0]=b[1]=5.0;
   } else if (icase <= 4) {
     Phia=0.0000001;
+    Phi[0]=0.0000001;
     Kab=Kac=Kad=-0.0004;
     La=1e-12;
     va=100.0;
     ba=5.0;
+    L[0]=1e-12;
+    v[0]=100.0;
+    b[0]=5.0;
+  }
 #else
   if (icase <= 1) {
     Na=Nb=0.0;
     Phia=Phib=0.0;
+    N[0]=N[1]=0.0;
+    Phi[0]=Phi[1]=0.0;
     Kab=Kac=Kad=Kbc=Kbd=0.0;
     La=Lb=Ld;
     va=vb=vd;
     ba=bb=0.0;
+    L[0]=L[1]=L[3];
+    v[0]=v[1]=v[3];
+    b[0]=b[1]=0.0;
   } else if (icase <= 4) {
     Na = 0.0;
     Phia=0.0;
+    N[0] = 0.0;
+    Phi[0]=0.0;
     Kab=Kac=Kad=0.0;
     La=Ld;
     va=vd;
     ba=0.0;
+    L[0]=L[3];
+    v[0]=v[3];
+    b[0]=0.0;
+  }
 #endif
   const double Xa = q*q*ba*ba*Na/6.0;
   const double Xb = q*q*bb*bb*Nb/6.0;
   const double Xc = q*q*bc*bc*Nc/6.0;
   const double Xd = q*q*bd*bd*Nd/6.0;
+  const double Xa = q*q*b[0]*b[0]*N[0]/6.0;
+  const double Xb = q*q*b[1]*b[1]*N[1]/6.0;
+  const double Xc = q*q*b[2]*b[2]*N[2]/6.0;
+  const double Xd = q*q*b[3]*b[3]*N[3]/6.0;
   // limit as Xa goes to 0 is 1
 …
 #if 0
   const double S0aa = icase<5
                       ? 1.0 : Na*Phia*va*Paa;
+                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa;
   const double S0bb = icase<2
                       ? 1.0 : Nb*Phib*vb*Pbb;
   const double S0cc = Nc*Phic*vc*Pcc;
   const double S0dd = Nd*Phid*vd*Pdd;
+                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb;
+  const double S0cc = N[2]*Phi[2]*v[2]*Pcc;
+  const double S0dd = N[3]*Phi[3]*v[3]*Pdd;
   const double S0ab = icase<8
                       ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
   const double S0ac = icase<9
                       ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
   const double S0ad = icase<9
                       ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
+                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
   const double S0bc = (icase!=4 && icase!=7 && icase!= 9)
                       ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
+                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
   const double S0bd = (icase!=4 && icase!=7 && icase!= 9)
                       ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
+                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
   const double S0cd = (icase==0 || icase==2 || icase==5)
                       ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
 #else  // sasview equivalent
 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc);
   double S0aa = Na*Phia*va*Paa;
   double S0bb = Nb*Phib*vb*Pbb;
   double S0cc = Nc*Phic*vc*Pcc;
   double S0dd = Nd*Phid*vd*Pdd;
   double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
   double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
   double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
   double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
   double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
   double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc);
+  double S0aa = N[0]*Phi[0]*v[0]*Paa;
+  double S0bb = N[1]*Phi[1]*v[1]*Pbb;
+  double S0cc = N[2]*Phi[2]*v[2]*Pcc;
+  double S0dd = N[3]*Phi[3]*v[3]*Pdd;
+  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
+  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
+  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
+  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
+  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
+  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
 switch(icase){
   case 0:
 …
   // Note: 1e-13 to convert from fm to cm for scattering length
   const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13;
   const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav;
   const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav;
   const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav;
+  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav;
+  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav;
+  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav;
   const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33
 …
+}
-double Iqxy(double qx, double qy,
-    double case_num,
-    double Na, double Phia, double va, double a_sld, double ba,
-    double Nb, double Phib, double vb, double b_sld, double bb,
-    double Nc, double Phic, double vc, double c_sld, double bc,
-    double Nd, double Phid, double vd, double d_sld, double bd,
-    double Kab, double Kac, double Kad,
-    double Kbc, double Kbd, double Kcd
+    )
+{
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q,
-        case_num,
-        Na, Phia, va, a_sld, ba,
-        Nb, Phib, vb, b_sld, bb,
-        Nc, Phic, vc, c_sld, bc,
-        Nd, Phid, vd, d_sld, bd,
-        Kab, Kac, Kad,
-        Kbc, Kbd, Kcd);
+}

sasmodels/models/rpa.py

-                      rec45c4f
+                      rea05c87
 #   ["name", "units", default, [lower, upper], "type","description"],
 parameters = [
     ["case_num", CASES, 0, [0, 10], "", "Component organization"],
+    ["case_num", "", 1, CASES, "", "Component organization"],
+    ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phia", "", 0.25, [0, 1], "", "volume fraction"],
+    ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["La", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["ba", "Ang", 5.0, [0, inf], "", "segment length"],
+    ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phib", "", 0.25, [0, 1], "", "volume fraction"],
+    ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["bb", "Ang", 5.0, [0, inf], "", "segment length"],
+    ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phic", "", 0.25, [0, 1], "", "volume fraction"],
+    ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["bc", "Ang", 5.0, [0, inf], "", "segment length"],
+    ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phid", "", 0.25, [0, 1], "", "volume fraction"],
+    ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["bd", "Ang", 5.0, [0, inf], "", "segment length"],
+    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"],
+    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"],
     ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"],

sasmodels/models/spherical_sld.py

-                      rec45c4f
+                      rd2bb604
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
 parameters = [["n_shells",         "",               1,      [0, 9],         "", "number of shells"],
+parameters = [["n",                "",               1,      [0, 9],         "", "number of shells"],
               ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"],
               ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"],
 …
 demo = dict(
     n_shells=4,
+    n=4,
     scale=1.0,
     solvent_sld=1.0,

sasmodels/models/squarewell.py

-                      rec45c4f
+                      rd2bb604
 """
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/stickyhardsphere.py

-                      rec45c4f
+                      rd2bb604
 """
-Iqxy = """
-    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS);
-    """
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/product.py

-                      rf247314
+                      rea05c87
 from .core import call_ER_VR
-from .generate import process_parameters
 SCALE=0
 …
     # Iq, Iqxy, form_volume, ER, VR and sesans
     model_info['composition'] = ('product', [p_info, s_info])
-    process_parameters(model_info)
     return model_info
 …
         # a parameter map.
         par_map = {}
         p_info = p_kernel.info['partype']
         s_info = s_kernel.info['partype']
+        p_info = p_kernel.info['par_type']
+        s_info = s_kernel.info['par_type']
         vol_pars = set(p_info['volume'])
         if dim == '2d':

sasmodels/resolution.py

-                      r4d76711
+                      rd2bb604
     from scipy.integrate import romberg
+    if any(k not in form.info['defaults'] for k in pars.keys()):
+        keys = set(form.info['defaults'].keys())
+        extra = set(pars.keys()) - keys
+        raise ValueError("bad parameters: [%s] not in [%s]"%
+                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))
+    par_set = set([p.name for p in form.info['parameters'].call_parameters])
+    if any(k not in par_set for k in pars.keys()):
+        extra = set(pars.keys()) - par_set
+        raise ValueError("bad parameters: [%s] not in [%s]"
+                         % (", ".join(sorted(extra)),
+                            ", ".join(sorted(pars.keys()))))
     if np.isscalar(width):
 …
     from scipy.integrate import romberg
+    if any(k not in form.info['defaults'] for k in pars.keys()):
+        keys = set(form.info['defaults'].keys())
+        extra = set(pars.keys()) - keys
+        raise ValueError("bad parameters: [%s] not in [%s]"%
+                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))
+    par_set = set([p.name for p in form.info['parameters'].call_parameters])
+    if any(k not in par_set for k in pars.keys()):
+        extra = set(pars.keys()) - par_set
+        raise ValueError("bad parameters: [%s] not in [%s]"
+                         % (", ".join(sorted(extra)),
+                            ", ".join(sorted(pars.keys()))))
     _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)

sasmodels/sasview_model.py

-                      r4d76711
+                      ree8f734
 from . import custom
 from . import generate
+from . import weights
 def load_standard_models():
 …
         try:
             models.append(_make_standard_model(name))
         except:
+        except Exception:
             logging.error(traceback.format_exc())
     return models
 …
         self._model = None
         model_info = self._model_info
+        parameters = model_info['parameters']
         self.name = model_info['name']
         self.description = model_info['description']
         self.category = None
+        self.multiplicity_info = None
+        self.is_multifunc = False
+        #self.is_multifunc = False
+        for p in parameters.kernel_parameters:
+            if p.is_control:
+                profile_axes = model_info['profile_axes']
+                self.multiplicity_info = [
+                    p.limits[1], p.name, p.choices, profile_axes[0]
+                    ]
+                break
+        else:
+            self.multiplicity_info = []
         ## interpret the parameters
 …
         self.params = collections.OrderedDict()
         self.dispersion = dict()
+        partype = model_info['partype']
+        for p in model_info['parameters']:
+        self.orientation_params = []
+        self.magnetic_params = []
+        self.fixed = []
+        for p in parameters.user_parameters():
             self.params[p.name] = p.default
             self.details[p.name] = [p.units] + p.limits
+        for name in partype['pd-2d']:
+            self.dispersion[name] = {
+                'width': 0,
+                'npts': 35,
+                'nsigmas': 3,
+                'type': 'gaussian',
+            }
+        self.orientation_params = (
+            partype['orientation']
+            + [n + '.width' for n in partype['orientation']]
+            + partype['magnetic'])
+        self.magnetic_params = partype['magnetic']
+        self.fixed = [n + '.width' for n in partype['pd-2d']]
+            if p.polydisperse:
+                self.dispersion[p.name] = {
+                    'width': 0,
+                    'npts': 35,
+                    'nsigmas': 3,
+                    'type': 'gaussian',
+                }
+            if p.type == 'orientation':
+                self.orientation_params.append(p.name)
+                self.orientation_params.append(p.name+".width")
+                self.fixed.append(p.name+".width")
+            if p.type == 'magnetic':
+                self.orientation_params.append(p.name)
+                self.magnetic_params.append(p.name)
+                self.fixed.append(p.name+".width")
         self.non_fittable = []
 …
         """
         # TODO: fix test so that parameter order doesn't matter
+        ret = ['%s.%s' % (d.lower(), p)
+               for d in self._model_info['partype']['pd-2d']
+               for p in ('npts', 'nsigmas', 'width')]
+        ret = ['%s.%s' % (p.name.lower(), ext)
+               for p in self._model_info['parameters'].user_parameters()
+               for ext in ('npts', 'nsigmas', 'width')
+               if p.polydisperse]
         #print(ret)
         return ret
 …
             # Check whether we have a list of ndarrays [qx,qy]
             qx, qy = qdist
+            partype = self._model_info['partype']
+            if not partype['orientation'] and not partype['magnetic']:
+            if not self._model_info['parameters'].has_2d:
                 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2))
             else:
 …
             self._model = core.build_model(self._model_info)
         q_vectors = [np.asarray(q) for q in args]
+        fn = self._model.make_kernel(q_vectors)
+        pars = [self.params[v] for v in fn.fixed_pars]
+        pd_pars = [self._get_weights(p) for p in fn.pd_pars]
+        result = fn(pars, pd_pars, self.cutoff)
+        fn.q_input.release()
+        fn.release()
+        kernel = self._model.make_kernel(q_vectors)
+        pairs = [self._get_weights(p)
+                 for p in self._model_info['parameters'].call_parameters]
+        details, weights, values = core.build_details(kernel, pairs)
+        result = kernel(details, weights, values, cutoff=self.cutoff)
+        kernel.q_input.release()
+        kernel.release()
         return result
 …
         Return dispersion weights for parameter
         """
+        from . import weights
+        relative = self._model_info['partype']['pd-rel']
+        limits = self._model_info['limits']
+        dis = self.dispersion[par]
+        value, weight = weights.get_weights(
+            dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
+            self.params[par], limits[par], par in relative)
+        return value, weight / np.sum(weight)
+        if par.polydisperse:
+            dis = self.dispersion[par.name]
+            value, weight = weights.get_weights(
+                dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
+                self.params[par.name], par.limits, par.relative_pd)
+            return value, weight / np.sum(weight)
+        else:
+            return [self.params[par.name]], []
 def test_model():

Changeset 3543141 in sasmodels

Legend:

Download in other formats: