Index: doc/developer/calculator.rst
===================================================================
--- doc/developer/calculator.rst	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
+++ doc/developer/calculator.rst	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
@@ -0,0 +1,201 @@
+Calculator Interface
+====================
+
+The environment needs to provide the following #defines:
+
+- USE_OPENCL is defined if running in opencl
+- KERNEL declares a function to be available externally
+- KERNEL_NAME is the name of the function being declared
+- NPARS is the number of parameters in the kernel
+- PARAMETER_TABLE is the declaration of the parameters to the kernel::
+
+    Cylinder:
+
+        #define PARAMETER_TABLE \
+        double length; \
+        double radius; \
+        double sld; \
+        double sld_solvent
+
+    Note: scale and background are never included
+
+    Multi-shell cylinder (10 shell max):
+
+        #define PARAMETER_TABLE \
+        double num_shells; \
+        double length; \
+        double radius[10]; \
+        double sld[10]; \
+        double sld_solvent
+
+- CALL_IQ(q, nq, i, pars) is the declaration of a call to the kernel::
+
+    Cylinder:
+
+        #define CALL_IQ(q, nq, i, var) \
+        Iq(q[i], \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent)
+
+    Multi-shell cylinder:
+
+        #define CALL_IQ(q, nq, i, var) \
+        Iq(q[i], \
+        var.num_shells, \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent)
+
+- CALL_VOLUME(var) is similar, but for calling the form volume::
+
+        #define CALL_VOLUME(var) \
+        form_volume(var.length, var.radius)
+
+- INVALID(var) is a test for model parameters in the correct range::
+
+    Cylinder:
+
+        #define INVALID(var) 0
+
+    BarBell:
+
+        #define INVALID(var) (var.bell_radius > var.radius)
+
+    Model with complicated constraints:
+
+        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
+we need to cycle through the values of the polydispersity, calculate
+the I(q, p) for each combination of parameters, and perform a normalized
+weighted sum across all the weights.  Parameters may be passed to the
+underlying calculation engine as scalars or vectors, but the polydispersity
+calculator treats the parameter set as one long vector.
+
+Let's assume we have 6 parameters in the model, with two polydisperse::
+
+    0: scale        {scl = constant}
+    1: background   {bkg = constant}
+    5: length       {l = vector of 30pts}
+    4: radius       {r = vector of 10pts}
+    3: sld          {s = constant/(radius**2*length)}
+    2: sld_solvent  {s2 = constant}
+
+This generates the following call to the kernel (where x stands for an
+arbitrary value that is not used by the kernel evaluator)::
+
+    NPARS = 4  // scale and background are in all models
+    problem {
+        pd_par = {5, 4, x, x}         // parameters *radius* and *length* vary
+        pd_length = {30, 10, 0, 0}    // *length* has more, so it is first
+        pd_offset = {10, 0, x, x}     // *length* starts at index 10 in weights
+        pd_stride = {1, 30, 300, 300} // cumulative product of pd length
+        pd_isvol = {1, 1, x, x}       // true if weight is a volume weight
+        par_offset = {2, 3, 303, 313}  // parameter offsets
+        par_coord = {0, 3, 2, 1} // bitmap of parameter dependencies
+        fast_coord_index = {5, 3, x, x}
+        fast_coord_count = 2  // two parameters vary with *length* distribution
+        theta_var = -1   // no spherical correction
+        fast_theta = 0   // spherical correction angle is not pd 1
+    }
+
+    weight = { l0, .., l29, r0, .., r9}
+    pars = { scl, bkg, l0, ..., l29, r0, r1, ..., r9,
+             s[l0,r0], ... s[l0,r9], s[l1,r0], ... s[l29,r9] , s2}
+
+    nq = 130
+    q = { q0, q1, ..., q130, x, x }  # pad to 8 element boundary
+    result = {r1, ..., r130, norm, vol, vol_norm, x, x, x, x, x, x, x}
+
+
+The polydisperse parameters are stored in as an array of parameter
+indices, one for each polydisperse parameter, stored in pd_par[n].
+Non-polydisperse parameters do not appear in this array. Each polydisperse
+parameter has a weight vector whose length is stored in pd_length[n],
+The weights are stored in a contiguous vector of weights for all
+parameters, with the starting position for the each parameter stored
+in pd_offset[n].  The values corresponding to the weights are stored
+together in a separate weights[] vector, with offset stored in
+par_offset[pd_par[n]]. Polydisperse parameters should be stored in
+decreasing order of length for highest efficiency.
+
+We limit the number of polydisperse dimensions to MAX_PD (currently 4).
+This cuts the size of the structure in half compared to allowing a
+separate polydispersity for each parameter.  This will help a little
+bit for models with large numbers of parameters, such as the onion model.
+
+Parameters may be coordinated.  That is, we may have the value of one
+parameter depend on a set of other parameters, some of which may be
+polydisperse.  For example, if sld is inversely proportional to the
+volume of a cylinder, and the length and radius are independently
+polydisperse, then for each combination of length and radius we need a
+separate value for the sld.  The caller must provide a coordination table
+for each parameter containing the value for each parameter given the
+value of the polydisperse parameters v1, v2, etc.  The tables for each
+parameter are arranged contiguously in a vector, with offset[k] giving the
+starting location of parameter k in the vector.  Each parameter defines
+coord[k] as a bit mask indicating which polydispersity parameters the
+parameter depends upon. Usually this is zero, indicating that the parameter
+is independent, but for the cylinder example given, the bits for the
+radius and length polydispersity parameters would both be set, the result
+being a (#radius x #length) table, or maybe a (#length x #radius) table
+if length comes first in the polydispersity table.
+
+NB: If we can guarantee that a compiler and OpenCL driver are available,
+we could instead create the coordination function on the fly for each
+parameter, saving memory and transfer time, but requiring a C compiler
+as part of the environment.
+
+In ordering the polydisperse parameters by decreasing length we can
+iterate over the longest dispersion weight vector first.  All parameters
+coordinated with this weight vector (the 'fast' parameters), can be
+updated with a simple increment to the next position in the parameter
+value table.  The indices of these parameters is stored in fast_coord_index[],
+with fast_coord_count being the number of fast parameters.  A total
+of NPARS slots is allocated to allow for the case that all parameters
+are coordinated with the fast index, though this will likely be mostly
+empty.  When the fast increment count reaches the end of the weight
+vector, then the index of the second polydisperse parameter must be
+incremented, and all of its coordinated parameters updated.  Because this
+operation is not in the inner loop, a slower algorithm can be used.
+
+If there is no polydispersity we pretend that it is polydisperisty with one
+parameter, pd_start=0 and pd_stop=1.  We may or may not short circuit the
+calculation in this case, depending on how much time it saves.
+
+The problem details structure can be allocated and sent in as an integer
+array using the read-only flag.  This allows us to copy it once per fit
+along with the weights vector, since features such as the number of
+polydisperity elements per pd parameter or the coordinated won't change
+between function evaluations.  A new parameter vector is sent for
+each I(q) evaluation.
+
+To protect against expensive evaluations taking all the GPU resource
+on large fits, the entire polydispersity will not be computed at once.
+Instead, a start and stop location will be sent, indicating where in the
+polydispersity loop the calculation should start and where it should
+stop.  We can do this for arbitrary start/stop points since we have
+unwound the nested loop.  Instead, we use the same technique as array
+index translation, using div and mod to figure out the i,j,k,...
+indices in the virtual nested loop.
+
+The results array will be initialized to zero for polydispersity loop
+entry zero, and preserved between calls to [start, stop] so that the
+results accumulate by the time the loop has completed.  Background and
+scale will be applied when the loop reaches the end.  This does require
+that the results array be allocated read-write, which is less efficient
+for the GPU, but it makes the calling sequence much more manageable.
+
+Scale and background cannot be coordinated with other polydisperse parameters
+
+Oriented objects in 2-D need a spherical correction on the angular variation
+in order to preserve the 'surface area' of the weight distribution.
+
+TODO: cutoff
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision 9ef9dd9f11c54d995dd94b650a4563230b868352)
+++ sasmodels/generate.py	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
@@ -21,4 +21,8 @@
 
     *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
@@ -216,5 +220,5 @@
 import sys
 from os.path import abspath, dirname, join as joinpath, exists, basename, \
-    splitext
+    splitext, getmtime
 import re
 import string
@@ -224,4 +228,5 @@
 import numpy as np
 
+# TODO: promote Parameter and model_info to classes
 PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description']
 Parameter = namedtuple('Parameter', PARAMETER_FIELDS)
@@ -230,5 +235,5 @@
 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
 
-C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
+TEMPLATE_ROOT = dirname(__file__)
 
 F16 = np.dtype('float16')
@@ -338,4 +343,5 @@
     raise ValueError("%r not found in %s" % (filename, search_path))
 
+
 def model_sources(model_info):
     """
@@ -346,17 +352,4 @@
     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 convert_type(source, dtype):
@@ -369,5 +362,5 @@
     if dtype == F16:
         fbytes = 2
-        source = _F16_PRAGMA + _convert_type(source, "half", "f")
+        source = _convert_type(source, "float", "f")
     elif dtype == F32:
         fbytes = 4
@@ -375,5 +368,5 @@
     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
@@ -418,28 +411,39 @@
 
 
-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())
+    return _template_cache[filename][1]
+
+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, ...) {
+             ....
+         }
+    """
+    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
+    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)
+
 def make_source(model_info):
     """
@@ -461,82 +465,47 @@
 
     # 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()
+    source = [load_template('kernel_header.c')]
 
     # Load additional sources
-    source = [open(f).read() for f in model_sources(model_info)]
+    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)]
+                     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 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 p.name in model_info['par_set']['2d']]
+    volume_parameters = model_info['par_type']['volume']
+
+    # Generate form_volume function, etc. from body only
+    if model_info['form_volume'] is not None:
+        pnames = [p.name for p in volume_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]
+        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]
+        source.append(_gen_fn('Iqxy', pnames, model_info['Iqxy']))
+
+    # Fill in definitions 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
-    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)
+        deref_vol = ",".join("v."+p.name for p in volume_parameters)
+        defines.append(('CALL_VOLUME(v)', 'form_volume(%s)\n'%deref_vol))
+    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(('IQ_KERNEL_NAME', model_info['name'] + '_Iq'))
+    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)))
-    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')))
-    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'))
@@ -557,15 +526,4 @@
         defines.append(('IQXY_CLOSE_LOOPS',
                         close_loops.replace('\n', ' \\\n')))
-    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
@@ -584,4 +542,29 @@
 def categorize_parameters(pars):
     """
+    Categorize the parameters by use:
+
+    * *pd* list of polydisperse parameters in order; gui should test whether
+      they are in *2d* or *magnetic* as appropriate for the data
+    * *1d* set of parameters that are used to compute 1D patterns
+    * *2d* set of parameters that are used to compute 2D patterns (which
+      includes all 1D parameters)
+    * *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
+      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']
+    return par_set
+
+def collect_types(pars):
+    """
     Build parameter categories out of the the parameter definitions.
 
@@ -596,38 +579,16 @@
     * *orientation* list of orientation parameters
     * *magnetic* list of magnetic parameters
-    * *<empty string>* list of parameters that have no type info
+    * *sld* list of parameters that have no type info
+    * *other* 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(),
+    """
+    par_type = {
+        'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], 'other': [],
     }
-
     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
+        par_type[p.type if p.type else 'other'].append(p.name)
+    return  par_type
+
 
 def process_parameters(model_info):
@@ -647,5 +608,6 @@
     partype = categorize_parameters(pars)
     model_info['limits'] = dict((p.name, p.limits) for p in pars)
-    model_info['partype'] = partype
+    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)
     if model_info.get('demo', None) is None:
Index: sasmodels/kernel_iq.c
===================================================================
--- sasmodels/kernel_iq.c	(revision 4b2972c83099d4d8433a1227a93599bde5b4422e)
+++ sasmodels/kernel_iq.c	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
@@ -182,5 +182,4 @@
       {
         const double scattering = CALL_IQ(q, nq, i, local_pars);
-        //const double scattering = Iq(q[i], IQ_PARAMETERS);
         result[i] += weight*scattering;
       }
@@ -200,5 +199,5 @@
         result[i] *= norm_vol/vol;
       }
-        result[i] = pars[0]*result[i]/norm + pars[1];
+      result[i] = pars[0]*result[i]/norm + pars[1];
     }
   }
Index: sasmodels/kernelpy.py
===================================================================
--- sasmodels/kernelpy.py	(revision 35647abd3c79d3ca6880d3b9913181dffce3a86f)
+++ sasmodels/kernelpy.py	(revision d5ac45fbf8e0b8649deccb5e4e72b31f7017e35b)
@@ -264,5 +264,5 @@
     result = scale * ret / norm + background
     return result
-=======
+
 """
 Python driver for python kernels
