Changes in / [31d22de:c2e4be1] in sasmodels

Files:

• doc/developer/calculator.rst (modified) (6 diffs)
• sasmodels/bumps_model.py (modified) (2 diffs)
• sasmodels/compare.py (modified) (2 diffs)
• sasmodels/core.py (modified) (7 diffs)
• sasmodels/direct_model.py (modified) (2 diffs)
• sasmodels/generate.py (modified) (12 diffs)
• sasmodels/kernel_header.c (modified) (4 diffs)
• sasmodels/kernel_iq.c (modified) (10 diffs)
• sasmodels/kerneldll.py (modified) (7 diffs)
• sasmodels/models/cylinder.c (modified) (3 diffs)
• sasmodels/models/gel_fit.c (modified) (4 diffs)
• sasmodels/product.py (modified) (1 diff)
• sasmodels/resolution.py (modified) (2 diffs)

doc/developer/calculator.rst

@@ -7 +7 @@
 - KERNEL declares a function to be available externally
 - KERNEL_NAME is the name of the function being declared
+- MAX_PD is the maximum depth of the polydispersity loop
 - NPARS is the number of parameters in the kernel
 - PARAMETER_TABLE is the declaration of the parameters to the kernel::
@@ -29 +30 @@
         double sld_solvent
 
-- CALL_IQ(q, nq, i, pars) is the declaration of a call to the kernel::
+- CALL_IQ(q, i, pars) is the declaration of a call to the kernel::
 
     Cylinder:
 
-        #define CALL_IQ(q, nq, i, var) \
-        Iq(q[i], \
+        #define CALL_IQ(q, i, var) Iq(q[i], \
         var.length, \
         var.radius, \
@@ -42 +42 @@
     Multi-shell cylinder:
 
-        #define CALL_IQ(q, nq, i, var) \
-        Iq(q[i], \
+        #define CALL_IQ(q, i, var) Iq(q[i], \
         var.num_shells, \
         var.length, \
@@ -50 +49 @@
         var.sld_solvent)
 
+    Cylinder2D:
+
+        #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent, \
+        var.theta, \
+        var.phi)
+
 - CALL_VOLUME(var) is similar, but for calling the form volume::
 
@@ -69 +78 @@
         inline bool constrained(p1, p2, p3) { return expression; }
         #define INVALID(var) constrained(var.p1, var.p2, var.p3)
-
-- IQ_FUNC could be Iq or Iqxy
-- IQ_PARS could be q[i] or q[2*i],q[2*i+1]
 
 Our design supports a limited number of polydispersity loops, wherein
@@ -200 +206 @@
 
 TODO: cutoff
+
+For accuracy we may want to introduce Kahan summation into the integration::
+
+
+    double accumulated_error = 0.0;
+    ...
+    #if USE_KAHAN_SUMMATION
+        const double y = next - accumulated_error;
+        const double t = ret + y;
+        accumulated_error = (t - ret) - y;
+        ret = t;
+    #else
+        ret += next;
+    #endif

sasmodels/bumps_model.py

@@ -81 +81 @@
     from bumps.names import Parameter
 
-    pars = {}
+    pars = {} # => floating point parameters
     for p in model_info['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 name in model_info['par_type']['pd']:
         for xpart, xdefault, xlimits in [
                 ('_pd', 0., pars[name].limits),
@@ -95 +95 @@
             pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits)
 
-    pd_types = {}
-    for name in model_info['partype']['pd-2d']:
+    pd_types = {}  # => distribution names
+    for name in model_info['par_type']['pd']:
         xname = name + '_pd_type'
         xvalue = kwargs.pop(xname, 'gaussian')

sasmodels/compare.py

@@ -308 +308 @@
         exclude = lambda n: False
     else:
-        partype = model_info['partype']
-        par1d = set(partype['fixed-1d']+partype['pd-1d'])
+        par1d = model_info['par_type']['1d']
         exclude = lambda n: n not in par1d
     lines = []
@@ -870 +869 @@
         pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars'])
         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 name in model_info['par_type']['1d']:
+                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
+                    k = 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/core.py

@@ -173 +173 @@
     return model(q_vectors)
 
-def get_weights(model_info, pars, name):
+def get_weights(parameter, values):
     """
     Generate the distribution for parameter *name* given the parameter values
@@ -181 +181 @@
     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 = values.get(parameter.name, parameter.default)
+    if parameter.type not in ('volume', 'orientation'):
+        return [value], []
+    relative = parameter.type == 'volume'
+    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)
     value, weight = weights.get_weights(
         disperser, npts, width, nsigma, value, limits, relative)
@@ -206 +208 @@
     return value, weight
 
-def call_kernel(kernel, pars, cutoff=0, mono=False):
+def call_kernel(kernel, values, cutoff=0, mono=False):
     """
     Call *kernel* returned from :func:`make_kernel` with parameters *pars*.
@@ -219 +221 @@
     *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]
-    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)
+    if mono or True:
+        pars = np.array([values.get(p.name, p.default)
+                         for p in kernel.info['parameters']])
+        weights = np.array([1.0])
+        details = kernel.info['mono_details']
+        return kernel(pars, weights, details, cutoff)
+    else:
+        pairs = [get_weights(p, values) for p in kernel.info['parameters']]
+        weights, pars = [v for v in zip(*pairs)]
+        details = generate.poly_details(kernel.info, weights, pars)
+        weights, pars = [np.hstack(v) for v in (weights, pars)]
+        return kernel(pars, weights, details, cutoff)
 
 def call_ER_VR(model_info, vol_pars):
@@ -249 +255 @@
 
 
-def call_ER(info, pars):
-    """
-    Call the model ER function using *pars*.
-    *info* is either *model.info* if you have a loaded model, or *kernel.info*
-    if you have a model kernel prepared for evaluation.
-    """
-    ER = info.get('ER', None)
+def call_ER(model_info, values):
+    """
+    Call the model ER function using *values*. *model_info* is either
+    *model.info* if you have a loaded model, or *kernel.info* if you
+    have a model kernel prepared for evaluation.
+    """
+    ER = model_info.get('ER', None)
     if ER is None:
         return 1.0
     else:
-        vol_pars = [get_weights(info, pars, name)
-                    for name in info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters']
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         individual_radii = ER(*value)
@@ -266 +273 @@
         return np.sum(weight*individual_radii) / np.sum(weight)
 
-def call_VR(info, pars):
+def call_VR(model_info, values):
     """
     Call the model VR function using *pars*.
@@ -272 +279 @@
     if you have a model kernel prepared for evaluation.
     """
-    VR = info.get('VR', None)
+    VR = model_info.get('VR', None)
     if VR is None:
         return 1.0
     else:
-        vol_pars = [get_weights(info, pars, name)
-                    for name in info['partype']['volume']]
+        vol_pars = [get_weights(parameter, values)
+                    for parameter in model_info['parameters']
+                    if parameter.type == 'volume']
         value, weight = dispersion_mesh(vol_pars)
         whole, part = VR(*value)

sasmodels/direct_model.py

@@ -66 +66 @@
             self.data_type = 'Iq'
 
-        partype = model.info['partype']
-
         if self.data_type == 'sesans':
             
@@ -81 +79 @@
             q_mono = sesans.make_all_q(data)
         elif self.data_type == 'Iqxy':
+            partype = model.info['par_type']
             if not partype['orientation'] and not partype['magnetic']:
                 raise ValueError("not 2D without orientation or magnetic parameters")

sasmodels/generate.py

@@ -236 +236 @@
 
 TEMPLATE_ROOT = dirname(__file__)
+
+MAX_PD = 4
 
 F16 = np.dtype('float16')
@@ -352 +354 @@
     return [_search(search_path, f) for f in model_info['source']]
 
-
 def convert_type(source, dtype):
     """
@@ -417 +418 @@
     if filename not in _template_cache or mtime > _template_cache[filename][0]:
         with open(path) as fid:
-            _template_cache[filename] = (mtime, fid.read())
+            _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 _gen_fn(name, pars, body):
@@ -431 +448 @@
          }
     """
-    template = """\
-double %(name)s(%(pars)s);
-double %(name)s(%(pars)s) {
-    %(body)s
-}
-
-
-"""
     par_decl = ', '.join('double ' + p for p in pars) if pars else 'void'
-    return template % {'name': name, 'body': body, 'pars': par_decl}
-
-def _gen_call_pars(name, pars):
-    name += "."
-    return ",".join(name+p for p in pars)
+    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl}
+
+def _call_pars(prefix, pars):
+    """
+    Return a list of *prefix.parameter* from parameter items.
+    """
+    prefix += "."
+    return [prefix+p 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):
@@ -464 +499 @@
     # for computing volume even if we allow non-disperse volume parameters.
 
-    # Load template
-    source = [load_template('kernel_header.c')]
-
-    # Load additional sources
-    source += [open(f).read() for f in model_sources(model_info)]
-
-    # Prepare defines
-    defines = []
-
-    iq_parameters = [p.name
-                     for p in model_info['parameters'][2:]  # skip scale, background
-                     if p.name in model_info['par_set']['1d']]
-    iqxy_parameters = [p.name
-                       for p in model_info['parameters'][2:]  # skip scale, background
-                       if p.name in model_info['par_set']['2d']]
-    volume_parameters = model_info['par_type']['volume']
+    # kernel_iq assumes scale and background are the first parameters;
+    # they should be first for 1d and 2d parameter lists as well.
+    assert model_info['parameters'][0].name == 'scale'
+    assert model_info['parameters'][1].name == 'background'
+
+    # Identify parameter types
+    iq_parameters = model_info['par_type']['1d'][2:]
+    iqxy_parameters = model_info['par_type']['2d'][2:]
+    vol_parameters = model_info['par_type']['volume']
+
+    # 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
 
     # Generate form_volume function, etc. from body only
     if model_info['form_volume'] is not None:
-        pnames = [p.name for p in volume_parameters]
+        pnames = [p for p in vol_parameters]
         source.append(_gen_fn('form_volume', pnames, model_info['form_volume']))
     if model_info['Iq'] is not None:
-        pnames = ['q'] + [p.name for p in iq_parameters]
+        pnames = ['q'] + [p for p in iq_parameters]
         source.append(_gen_fn('Iq', pnames, model_info['Iq']))
     if model_info['Iqxy'] is not None:
-        pnames = ['qx', 'qy'] + [p.name for p in iqxy_parameters]
+        pnames = ['qx', 'qy'] + [p for p in iqxy_parameters]
         source.append(_gen_fn('Iqxy', pnames, model_info['Iqxy']))
 
-    # Fill in definitions for volume parameters
-    if volume_parameters:
-        deref_vol = ",".join("v."+p.name for p in volume_parameters)
-        defines.append(('CALL_VOLUME(v)', 'form_volume(%s)\n'%deref_vol))
+    # Define the parameter table
+    source.append("#define PARAMETER_TABLE \\")
+    source.append("\\\n    ".join("double %s;"%p.name
+                                   for p in model_info['parameters'][2:]))
+
+    # Define the function calls
+    if vol_parameters:
+        refs = ",".join(_call_pars("v", vol_parameters))
+        call_volume = "#define CALL_VOLUME(v) form_volume(%s)"%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.
-        defines.append(('CALL_VOLUME(v)', '0.0'))
-
-    # Fill in definitions for Iq parameters
-    defines.append(('KERNEL_NAME', model_info['name']))
-    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)))
-    # 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')))
-    # 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,
-        }
+        call_volume = "#define CALL_VOLUME(v) 0.0"
+    source.append(call_volume)
+
+    refs = ["q[i]"] + _call_pars("v", 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", 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"%model_info['max_pd'])
+    source.append("#define NPARS %d"%(len(model_info['parameters'])-2))
+
+    # 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)
 
 def categorize_parameters(pars):
@@ -551 +590 @@
     * *magnetic* set of parameters that are used to compute magnetic
       patterns (which includes all 1D and 2D parameters)
-    * *sesans* set of parameters that are used to compute sesans patterns
-     (which is just 1D without background)
-    * *pd-relative* is the set of parameters with relative distribution
+    * *pd_relative* is the set of parameters with relative distribution
       width (e.g., radius +/- 10%) rather than absolute distribution
       width (e.g., theta +/- 6 degrees).
     """
     par_set = {}
-    par_set['1d'] = [p for p in pars if p.type not in ('orientation', 'magnetic')]
-    par_set['2d'] = [p for p in pars if p.type != 'magnetic']
-    par_set['magnetic'] = [p for p in pars]
-    par_set['pd'] = [p for p in pars if p.type in ('volume', 'orientation')]
-    par_set['pd_relative'] = [p for p in pars if p.type == 'volume']
+    par_set['1d'] = [p.name for p in pars if p.type not in ('orientation', 'magnetic')]
+    par_set['2d'] = [p.name for p in pars if p.type != 'magnetic']
+    par_set['magnetic'] = [p.name for p in pars]
+    par_set['pd'] = [p.name for p in pars if p.type in ('volume', 'orientation')]
+    par_set['pd_relative'] = [p.name for p in pars if p.type == 'volume']
     return par_set
 
@@ -605 +642 @@
     pars = [Parameter(*p) for p in model_info['parameters']]
     # Fill in the derived attributes
+    par_type = collect_types(pars)
+    par_type.update(categorize_parameters(pars))
     model_info['parameters'] = pars
-    partype = categorize_parameters(pars)
-    model_info['limits'] = dict((p.name, p.limits) for p in pars)
-    model_info['par_type'] = collect_types(pars)
-    model_info['par_set'] = categorize_parameters(pars)
-    model_info['defaults'] = dict((p.name, p.default) for p in pars)
+    model_info['par_type'] = par_type
     if model_info.get('demo', None) is None:
-        model_info['demo'] = model_info['defaults']
-    model_info['has_2d'] = partype['orientation'] or partype['magnetic']
+        model_info['demo'] = dict((p.name, p.default) for p in pars)
+    model_info['has_2d'] = par_type['orientation'] or par_type['magnetic']
+
+def mono_details(max_pd, npars):
+    par_offset = 5*max_pd
+    const_offset = par_offset + 3*npars
+
+    mono = np.zeros(const_offset + 3, 'i')
+    mono[0] = 0                # pd_par: arbitrary order; use first
+    mono[1*max_pd] = 1         # pd_length: only one element
+    mono[2*max_pd] = 2         # pd_offset: skip scale and background
+    mono[3*max_pd] = 1         # pd_stride: vectors of length 1
+    mono[4*max_pd-1] = 1       # pd_stride[-1]: only one element in set
+    mono[4*max_pd] = 0         # pd_isvol: doens't matter if no norm
+    mono[par_offset:par_offset+npars] = np.arange(2, npars+2, dtype='i')
+    # par_offset: copied in order
+    mono[par_offset+npars:par_offset+2*npars] = 0
+    # par_coord: no coordinated parameters
+    mono[par_offset+npars] = 1 # par_coord[0]: except par 0
+    mono[par_offset+2*npars:par_offset+3*npars] = 0
+    # fast coord with 0
+    mono[const_offset] = 1     # fast_coord_count: one fast index
+    mono[const_offset+1] = -1  # theta_var: None
+    mono[const_offset+2] = 0   # fast_theta: False
+    return mono
+
+def poly_details(model_info, weights, pars, constraints=None):
+    if constraints is not None:
+        # 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 ValueError("Can't handle constraints yet")
+
+    raise ValueError("polydispersity not supported")
+
+
+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.
+    """
+    if model_info['Iq'] is not None and model_info['Iqxy'] is None:
+        if model_info['par_type']['1d'] != model_info['par_type']['2d']:
+            raise ValueError("Iqxy model is missing")
+        Iq = model_info['Iq']
+        def Iqxy(qx, qy, **kw):
+            return Iq(np.sqrt(qx**2 + qy**2), **kw)
+        model_info['Iqxy'] = Iqxy
 
 def make_model_info(kernel_module):
@@ -640 +724 @@
       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
@@ -661 +741 @@
       *model_info* blocks for the composition objects.  This allows us to
       build complete product and mixture models from just the info.
+    * *max_pd* is the max polydispersity dimension.  This is constant and
+      should not be reset.  You may be able to change it when the program
+      starts by setting *sasmodels.generate.MAX_PD*.
 
     """
@@ -671 +754 @@
         name = " ".join(w.capitalize() for w in kernel_id.split('_'))
     model_info = dict(
+        max_pd=MAX_PD,
         id=kernel_id,  # string used to load the kernel
         filename=abspath(kernel_module.__file__),
@@ -688 +772 @@
         oldpars=getattr(kernel_module, 'oldpars', {}),
         tests=getattr(kernel_module, 'tests', []),
+        mono_details = mono_details(MAX_PD, len(kernel_module.parameters))
         )
     process_parameters(model_info)
@@ -693 +778 @@
     functions = "ER VR form_volume Iq Iqxy shape sesans".split()
     model_info.update((k, getattr(kernel_module, k, None)) for k in functions)
+    create_default_functions(model_info)
     return model_info
 

sasmodels/kernel_header.c

@@ -1 +1 @@
 #ifdef __OPENCL_VERSION__
 # define USE_OPENCL
+#elif defined(_OPENMP)
+# define USE_OPENMP
 #endif
-
-#define USE_KAHAN_SUMMATION 0
 
 // If opencl is not available, then we are compiling a C function
@@ -16 +16 @@
          #include <float.h>
          #define kernel extern "C" __declspec( dllexport )
-         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); }
-         inline double fmin(double x, double y) { return x>y ? y : x; }
-         inline double fmax(double x, double y) { return x<y ? y : x; }
-         inline double isnan(double x) { return _isnan(x); }
+         static double trunc(double x) { return x>=0?floor(x):-floor(-x); }
+         static double fmin(double x, double y) { return x>y ? y : x; }
+         static double fmax(double x, double y) { return x<y ? y : x; }
+         static double isnan(double x) { return _isnan(x); }
          #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
          static double cephes_expm1(double x) {
@@ -48 +48 @@
          #define kernel extern "C"
      #endif
-     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
+     static void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
 #  else
      #include <stdio.h>
@@ -99 +99 @@
 #  define M_4PI_3 4.18879020478639
 #endif
-//inline double square(double x) { return pow(x,2.0); }
-//inline double square(double x) { return pown(x,2); }
-inline double square(double x) { return x*x; }
-inline double cube(double x) { return x*x*x; }
-inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
+static double square(double x) { return x*x; }
+static double cube(double x) { return x*x*x; }
+static double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
 

sasmodels/kernel_iq.c

@@ -12 +12 @@
 */
 
+#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
+#define _PAR_BLOCK_
 
 #define MAX_PD 4  // MAX_PD is the max number of polydisperse parameters
-#define PD_2N 16  // PD_2N is the size of the coordination step table
 
 typedef struct {
@@ -31 +32 @@
 
 typedef struct {
-    PARAMETER_DECL;
+    PARAMETER_TABLE;
 } ParameterBlock;
-
-#define FULL_KERNEL_NAME KERNEL_NAME ## _ ## IQ_FUNC
-KERNEL
-void FULL_KERNEL_NAME(
+#endif
+
+
+kernel
+void KERNEL_NAME(
     int nq,                 // number of q values
+    const int pd_start,     // where we are in the polydispersity loop
+    const int pd_stop,      // where we are stopping in the polydispersity loop
     global const ProblemDetails *problem,
     global const double *weights,
     global const double *pars,
     global const double *q, // nq q values, with padding to boundary
-    global double *result,  // nq return values, again with padding
-    const double cutoff,    // cutoff in the polydispersity weight product
-    const int pd_start,     // where we are in the polydispersity loop
-    const int pd_stop,      // where we are stopping in the polydispersity loop
+    global double *result,  // nq+3 return values, again with padding
+    const double cutoff     // cutoff in the polydispersity weight product
     )
 {
-
+printf("hello\n");
   // Storage for the current parameter values.  These will be updated as we
   // walk the polydispersity cube.
   local ParameterBlock local_pars;  // current parameter values
-  const double *parvec = &local_pars;  // Alias named parameters with a vector
+  double *pvec = (double *)(&local_pars);  // Alias named parameters with a vector
 
   local int offset[NPARS-2];
@@ -60 +62 @@
     // Shouldn't need to copy!!
     for (int k=0; k < NPARS; k++) {
-      parvec[k] = pars[k+2];  // skip scale and background
+      pvec[k] = pars[k+2];  // skip scale and background
     }
 
@@ -68 +70 @@
     for (int i=0; i < nq; i++) {
     {
-      const double scattering = IQ_FUNC(IQ_PARS, IQ_PARAMETERS);
+      const double scattering = CALL_IQ(q, i, local_pars);
       result[i] += pars[0]*scattering + pars[1];
     }
@@ -99 +101 @@
     norm = result[nq];
     vol = result[nq+1];
-    norm_vol = results[nq+2];
+    norm_vol = result[nq+2];
   }
 
   // Location in the polydispersity hypercube, one index per dimension.
-  local int pd_index[PD_MAX];
+  local int pd_index[MAX_PD];
 
   // Trigger the reset behaviour that happens at the end the fast loop
   // by setting the initial index >= weight vector length.
-  pd_index[0] = pd_length[0];
+  pd_index[0] = problem->pd_length[0];
+
 
   // need product of weights at every Iq calc, so keep product of
@@ -120 +123 @@
   for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
     // check if indices need to be updated
-    if (pd_index[0] >= pd_length[0]) {
+    if (pd_index[0] >= problem->pd_length[0]) {
 
       // RESET INDICES
-      pd_index[0] = loop_index%pd_length[0];
+      pd_index[0] = loop_index%problem->pd_length[0];
       partial_weight = 1.0;
       partial_volweight = 1.0;
       for (int k=1; k < MAX_PD; k++) {
-        pd_index[k] = (loop_index%pd_length[k])/pd_stride[k];
-        const double wi = weights[pd_offset[0]+pd_index[0]];
+        pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k];
+        const double wi = weights[problem->pd_offset[0]+pd_index[0]];
         partial_weight *= wi;
-        if (pd_isvol[k]) partial_volweight *= wi;
+        if (problem->pd_isvol[k]) partial_volweight *= wi;
       }
       for (int k=0; k < NPARS; k++) {
-        int coord = par_coord[k];
-        int this_offset = par_offset[k];
+        int coord = problem->par_coord[k];
+        int this_offset = problem->par_offset[k];
         int block_size = 1;
         for (int bit=0; bit < MAX_PD && coord != 0; bit++) {
           if (coord&1) {
               this_offset += block_size * pd_index[bit];
-              block_size *= pd_length[bit];
+              block_size *= problem->pd_length[bit];
           }
           coord /= 2;
         }
         offset[k] = this_offset;
-        parvec[k] = pars[this_offset];
-      }
-      weight = partial_weight * weights[pd_offset[0]+pd_index[0]]
-      if (theta_var >= 0) {
-        spherical_correction = fabs(cos(M_PI_180*parvec[theta_var]));
-      }
-      if (!fast_theta) weight *= spherical_correction;
+        pvec[k] = pars[this_offset];
+      }
+      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]];
+      if (problem->theta_var >= 0) {
+        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_var]));
+      }
+      if (!problem->fast_theta) {
+        weight *= spherical_correction;
+      }
 
     } else {
@@ -156 +161 @@
       // INCREMENT INDICES
       pd_index[0] += 1;
-      const double wi = weights[pd_offset[0]+pd_index[0]];
+      const double wi = weights[problem->pd_offset[0]+pd_index[0]];
       weight = partial_weight*wi;
-      if (pd_isvol[0]) vol_weight *= wi;
+      if (problem->pd_isvol[0]) vol_weight *= wi;
       for (int k=0; k < problem->fast_coord_count; k++) {
-        parvec[ fast_coord_index[k]]
-            = pars[offset[fast_coord_index[k]] + pd_index[0]];
-      }
-      if (fast_theta) weight *= fabs(cos(M_PI_180*parvec[theta_var]));
-
+        pvec[problem->fast_coord_index[k]]
+            = pars[offset[problem->fast_coord_index[k]]++];
+      }
+      if (problem->fast_theta) {
+        weight *= fabs(cos(M_PI_180*pvec[problem->theta_var]));
+      }
     }
 
@@ -171 +177 @@
     #endif
 
+    // Accumulate I(q)
+    // Note: weight==0 must always be excluded
     if (weight > cutoff) {
       norm += weight;
@@ -180 +188 @@
       #endif
       for (int i=0; i < nq; i++) {
-      {
-        const double scattering = CALL_IQ(q, nq, i, local_pars);
+        const double scattering = CALL_IQ(q, i, local_pars);
         result[i] += weight*scattering;
       }
+    }
   }
   //Makes a normalization avialable for the next round
@@ -191 +199 @@
 
   //End of the PD loop we can normalize
-  if (pd_stop == pd_stride[MAX_PD-1]) {}
+  if (pd_stop >= problem->pd_stride[MAX_PD-1]) {
     #ifdef USE_OPENMP
     #pragma omp parallel for

sasmodels/kerneldll.py

@@ -61 +61 @@
 if sys.platform == 'darwin':
     #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
-    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
+    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm -fno-unknown-pragmas"
 elif os.name == 'nt':
     # call vcvarsall.bat before compiling to set path, headers, libs, etc.
@@ -80 +80 @@
         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 += " -fWno-unknown-pragmas"
 else:
     COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
@@ -140 +142 @@
 
     source = generate.convert_type(source, dtype)
-    source_files = generate.model_sources(model_info) + [model_info['filename']]
+    source_files = (generate.model_sources(model_info)
+                    + [model_info['filename']]
+                    + generate.model_templates())
     dll = dll_path(model_info, dtype)
     newest = max(os.path.getmtime(f) for f in source_files)
@@ -170 +174 @@
     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):
     """
@@ -195 +195 @@
 
     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:
@@ -210 +205 @@
               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_int]*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.Iq.argtypes = argtypes
+        self.Iqxy.argtypes = argtypes
 
     def __getstate__(self):
@@ -261 +256 @@
         self.q_input = q_input
         self.kernel = kernel
-        self.res = np.empty(q_input.nq, q_input.dtype)
+        self.res = np.empty(q_input.nq+3, 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.parameters = model_info['par_type'][dim]
 
         # In dll kernel, but not in opencl kernel
         self.p_res = self.res.ctypes.data
 
-    def __call__(self, fixed_pars, pd_pars, cutoff):
+    def __call__(self, details, values, weights, 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)
+        args = [
+            self.q_input.nq, # nq
+            0, # pd_start
+            1, # pd_stop
+            details.ctypes.data, # problem
+            weights.ctypes.data,  # weights
+            values.ctypes.data,  #pars
+            self.q_input.q_pointers[0], #q
+            self.p_res,   # results
+            real(cutoff), # cutoff
+            ]
         self.kernel(*args)
 
-        return self.res
+        return self.res[:-3]
 
     def release(self):

sasmodels/models/cylinder.c

@@ -3 +3 @@
 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)
@@ -15 +17 @@
     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;
@@ -47 +48 @@
     double phi)
 {
-    // TODO: return NaN if radius<0 or length<0?
     double sn, cn; // slots to hold sincos function output
 

sasmodels/models/gel_fit.c

@@ -1 @@
-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,
@@ -24 +8 @@
     // 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) );
@@ -30 +14 @@
     // 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));
 
@@ -37 +21 @@
     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/product.py

@@ -99 +99 @@
         # 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

@@ -504 +504 @@
     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
+    par_set = set([p.name for p in form.info['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(keys))))
@@ -558 +558 @@
     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
+    par_set = set([p.name for p in form.info['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(keys))))