Index: doc/developer/calculator.rst
===================================================================
--- doc/developer/calculator.rst	(revision 6fd8de01432581dc71f1290cf5070acac4860fd7)
+++ doc/developer/calculator.rst	(revision 6fd8de01432581dc71f1290cf5070acac4860fd7)
@@ -0,0 +1,225 @@
+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
+- MAX_PD is the maximum depth of the polydispersity loop [model specific]
+- 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, i, var) is the declaration of a call to the kernel::
+
+    Cylinder:
+
+        #define CALL_IQ(q, i, var) Iq(q[i], \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        var.sld_solvent)
+
+    Multi-shell cylinder:
+
+        #define CALL_IQ(q, i, var) Iq(q[i], \
+        var.num_shells, \
+        var.length, \
+        var.radius, \
+        var.sld, \
+        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::
+
+        #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)
+
+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}
+    2: length       {l = vector of 30pts}
+    3: radius       {r = vector of 10pts}
+    4: sld          {s = constant/(radius**2*length)}
+    5: 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 = {3, 2, 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 = {True, True, 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 = {3, 5, x, x} // radius and sld have fast index
+        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} //length comes first as the longest vec
+    pars = { scl, bkg, l0, ..., l29, r0, r1, ..., r9,
+             s[l0,r0], ... s[l0,r9], s[l1,r0], ... s[l29,r9] , s2}
+             //where s[x,y] stands for material sld, s2 = solvent sld
+
+    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),
+though some models may have fewer if they have fewer polydisperse
+parameters. 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
+par_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.
+
+cutoff parameter limits integration area within polydispersity hypercude,
+which speeds calculations
+
+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
Index: doc/developer/index.rst
===================================================================
--- doc/developer/index.rst	(revision 56fc97aa819cfe686101e6c6b3bac4c2c7cfa01e)
+++ doc/developer/index.rst	(revision b85be2da3b04a7ce8ef87e2a85a44572b2b1480c)
@@ -3,3 +3,7 @@
 ***************************
 
+.. toctree::
+   :numbered: 4
+   :maxdepth: 4
 
+   calculator.rst
Index: doc/genapi.py
===================================================================
--- doc/genapi.py	(revision 3d5c6f859c4a9840b418a6d68e9b567185d00b12)
+++ doc/genapi.py	(revision a5b8477c53695e7df03ecec7346f0c1bd770d3f7)
@@ -59,19 +59,28 @@
     #('alignment', 'GPU data alignment [unused]'),
     ('bumps_model', 'Bumps interface'),
+    ('compare', 'Compare models on different compute engines'),
     ('convert', 'Sasview to sasmodel converter'),
     ('core', 'Model access'),
+    ('data', 'Data layout and plotting routines'),
+    ('details', 'Parameter packing for kernel calls'),
     ('direct_model', 'Simple interface'),
     ('exception', 'Annotate exceptions'),
+    #('frozendict', 'Freeze a dictionary to make it immutable'),
     ('generate', 'Model parser'),
+    ('kernel', 'Evaluator type definitions'),
     ('kernelcl', 'OpenCL model evaluator'),
     ('kerneldll', 'Ctypes model evaluator'),
     ('kernelpy', 'Python model evaluator'),
+    ('list_pars', 'Identify all parameters in all models'),
+    ('mixture', 'Mixture model evaluator'),
     ('model_test', 'Unit test support'),
+    ('modelinfo', 'Parameter and model definitions'),
+    ('product', 'Product model evaluator'),
     ('resolution', '1-D resolution functions'),
     ('resolution2d', '2-D resolution functions'),
     ('sasview_model', 'Sasview interface'),
-    ('sesans', 'SESANS model evaluator'),
+    ('sesans', 'SESANS calculation routines'),
+    #('transition', 'Model stepper for automatic model selection'),
     ('weights', 'Distribution functions'),
-    #('transition', 'Model stepper for automatic model selection'),
 ]
 package='sasmodels'
Index: doc/genmodel.py
===================================================================
--- doc/genmodel.py	(revision 89f41635a25c3421f84f194bd2384d939489b0de)
+++ doc/genmodel.py	(revision a5b8477c53695e7df03ecec7346f0c1bd770d3f7)
@@ -1,6 +1,7 @@
+from __future__ import print_function
+
 import sys, os, math, re
 import numpy as np
 import matplotlib.pyplot as plt
-import pylab
 sys.path.insert(0, os.path.abspath('..'))
 from sasmodels import generate, core
@@ -8,35 +9,18 @@
 from sasmodels.data import empty_data1D, empty_data2D
 
-
-# Convert ../sasmodels/models/name.py to name
-model_name = os.path.basename(sys.argv[1])[:-3]
-model_info = core.load_model_info(model_name)
-model = core.build_model(model_info)
-
-# Load the doc string from the module definition file and store it in rst
-docstr = generate.make_doc(model_info)
-
-
-# Calculate 1D curve for default parameters
-pars = dict((p.name, p.default) for p in model_info['parameters'])
-
-# Plotting ranges and options
-opts = {
-    'xscale'    : 'log',
-    'yscale'    : 'log' if not model_info['structure_factor'] else 'linear',
-    'zscale'    : 'log' if not model_info['structure_factor'] else 'linear',
-    'q_min'     : 0.001,
-    'q_max'     : 1.0,
-    'nq'        : 1000,
-    'nq2d'      : 1000,
-    'vmin'      : 1e-3,  # floor for the 2D data results
-    'qx_max'    : 0.5,
-    #'colormap'  : 'gist_ncar',
-    'colormap'  : 'nipy_spectral',
-    #'colormap'  : 'jet',
-}
-
+try:
+    from typing import Dict, Any
+except ImportError:
+    pass
+else:
+    from matplotlib.axes import Axes
+    from sasmodels.kernel import KernelModel
+    from sasmodels.modelinfo import ModelInfo
 
 def plot_1d(model, opts, ax):
+    # type: (KernelModel, Dict[str, Any], Axes) -> None
+    """
+    Create a 1-D image.
+    """
     q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq']
     q_min = math.log10(q_min)
@@ -47,5 +31,5 @@
     Iq1D = calculator()
 
-    ax.plot(q, Iq1D, color='blue', lw=2, label=model_info['name'])
+    ax.plot(q, Iq1D, color='blue', lw=2, label=model.info.name)
     ax.set_xlabel(r'$Q \/(\AA^{-1})$')
     ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$')
@@ -55,6 +39,10 @@
 
 def plot_2d(model, opts, ax):
+    # type: (KernelModel, Dict[str, Any], Axes) -> None
+    """
+    Create a 2-D image.
+    """
     qx_max, nq2d = opts['qx_max'], opts['nq2d']
-    q = np.linspace(-qx_max, qx_max, nq2d)
+    q = np.linspace(-qx_max, qx_max, nq2d) # type: np.ndarray
     data2d = empty_data2D(q, resolution=0.0)
     calculator = DirectModel(data2d, model)
@@ -64,72 +52,123 @@
         Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf))
     ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='lower',
-        extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=opts['colormap'])
+              extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=opts['colormap'])
     ax.set_xlabel(r'$Q_x \/(\AA^{-1})$')
     ax.set_ylabel(r'$Q_y \/(\AA^{-1})$')
 
-# Generate image 
-fig_height = 3.0 # in
-fig_left = 0.6 # in
-fig_right = 0.5 # in
-fig_top = 0.6*0.25 # in
-fig_bottom = 0.6*0.75
-if model_info['has_2d']:
-    plot_height = fig_height - (fig_top+fig_bottom)
-    plot_width = plot_height
-    fig_width = 2*(plot_width + fig_left + fig_right)
-    aspect = (fig_width, fig_height)
-    ratio = aspect[0]/aspect[1]
-    ax_left = fig_left/fig_width
-    ax_bottom = fig_bottom/fig_height
-    ax_height = plot_height/fig_height
-    ax_width = ax_height/ratio # square axes
-    fig = plt.figure(figsize=aspect)
-    ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height])
-    plot_2d(model, opts, ax2d)
-    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
-    plot_1d(model, opts, ax1d)
-    #ax.set_aspect('square')
-else:
-    plot_height = fig_height - (fig_top+fig_bottom)
-    plot_width = (1+np.sqrt(5))/2*fig_height
-    fig_width = plot_width + fig_left + fig_right
-    ax_left = fig_left/fig_width
-    ax_bottom = fig_bottom/fig_height
-    ax_width = plot_width/fig_width
-    ax_height = plot_height/fig_height
-    aspect = (fig_width, fig_height)
-    fig = plt.figure(figsize=aspect)
-    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
-    plot_1d(model, opts, ax1d)
+def figfile(model_info):
+    # type: (ModelInfo) -> str
+    return model_info.id + '_autogenfig.png'
 
-# Save image in model/img
-figname = model_name + '_autogenfig.png'
-filename = os.path.join('model', 'img', figname)
-plt.savefig(filename, bbox_inches='tight')
-#print "figure saved in",filename
+def make_figure(model_info, opts):
+    # type: (ModelInfo, Dict[str, Any]) -> None
+    """
+    Generate the figure file to include in the docs.
+    """
+    model = core.build_model(model_info)
 
-# Auto caption for figure
-captionstr = '\n'
-captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n'
-captionstr += '\n'
-if model_info['has_2d']:
-    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n'
-else:
-    captionstr += '    1D plot corresponding to the default parameters of the model.\n'
-captionstr += '\n'
+    fig_height = 3.0 # in
+    fig_left = 0.6 # in
+    fig_right = 0.5 # in
+    fig_top = 0.6*0.25 # in
+    fig_bottom = 0.6*0.75
+    if model_info.parameters.has_2d:
+        plot_height = fig_height - (fig_top+fig_bottom)
+        plot_width = plot_height
+        fig_width = 2*(plot_width + fig_left + fig_right)
+        aspect = (fig_width, fig_height)
+        ratio = aspect[0]/aspect[1]
+        ax_left = fig_left/fig_width
+        ax_bottom = fig_bottom/fig_height
+        ax_height = plot_height/fig_height
+        ax_width = ax_height/ratio # square axes
+        fig = plt.figure(figsize=aspect)
+        ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height])
+        plot_2d(model, opts, ax2d)
+        ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
+        plot_1d(model, opts, ax1d)
+        #ax.set_aspect('square')
+    else:
+        plot_height = fig_height - (fig_top+fig_bottom)
+        plot_width = (1+np.sqrt(5))/2*fig_height
+        fig_width = plot_width + fig_left + fig_right
+        ax_left = fig_left/fig_width
+        ax_bottom = fig_bottom/fig_height
+        ax_width = plot_width/fig_width
+        ax_height = plot_height/fig_height
+        aspect = (fig_width, fig_height)
+        fig = plt.figure(figsize=aspect)
+        ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height])
+        plot_1d(model, opts, ax1d)
 
-# Add figure reference and caption to documentation (at end, before References)
-pattern = '\*\*REFERENCE'
-m = re.search(pattern, docstr.upper())
+    # Save image in model/img
+    path = os.path.join('model', 'img', figfile(model_info))
+    plt.savefig(path, bbox_inches='tight')
+    #print("figure saved in",path)
 
-if m:
-    docstr1 = docstr[:m.start()]
-    docstr2 = docstr[m.start():]
-    docstr = docstr1 + captionstr + docstr2
-else:
-    print '------------------------------------------------------------------'
-    print 'References NOT FOUND for model: ', model_info['id']
-    print '------------------------------------------------------------------'
-    docstr = docstr + captionstr
+def gen_docs(model_info):
+    # type: (ModelInfo) -> None
+    """
+    Generate the doc string with the figure inserted before the references.
+    """
 
-open(sys.argv[2],'w').write(docstr)
+    # Load the doc string from the module definition file and store it in rst
+    docstr = generate.make_doc(model_info)
+
+    # Auto caption for figure
+    captionstr = '\n'
+    captionstr += '.. figure:: img/' + figfile(model_info) + '\n'
+    captionstr += '\n'
+    if model_info.parameters.has_2d:
+        captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n'
+    else:
+        captionstr += '    1D plot corresponding to the default parameters of the model.\n'
+    captionstr += '\n'
+
+    # Add figure reference and caption to documentation (at end, before References)
+    pattern = '\*\*REFERENCE'
+    match = re.search(pattern, docstr.upper())
+
+    if match:
+        docstr1 = docstr[:match.start()]
+        docstr2 = docstr[match.start():]
+        docstr = docstr1 + captionstr + docstr2
+    else:
+        print('------------------------------------------------------------------')
+        print('References NOT FOUND for model: ', model_info.id)
+        print('------------------------------------------------------------------')
+        docstr += captionstr
+
+    open(sys.argv[2],'w').write(docstr)
+
+def process_model(path):
+    # type: (str) -> None
+    """
+    Generate doc file and image file for the given model definition file.
+    """
+
+    # Load the model file
+    model_name = os.path.basename(path)[:-3]
+    model_info = core.load_model_info(model_name)
+
+    # Plotting ranges and options
+    opts = {
+        'xscale'    : 'log',
+        'yscale'    : 'log' if not model_info.structure_factor else 'linear',
+        'zscale'    : 'log' if not model_info.structure_factor else 'linear',
+        'q_min'     : 0.001,
+        'q_max'     : 1.0,
+        'nq'        : 1000,
+        'nq2d'      : 1000,
+        'vmin'      : 1e-3,  # floor for the 2D data results
+        'qx_max'    : 0.5,
+        #'colormap'  : 'gist_ncar',
+        'colormap'  : 'nipy_spectral',
+        #'colormap'  : 'jet',
+    }
+
+    # Generate the RST file and the figure.  Order doesn't matter.
+    gen_docs(model_info)
+    make_figure(model_info, opts)
+
+if __name__ == "__main__":
+    process_model(sys.argv[1])
Index: doc/gentoc.py
===================================================================
--- doc/gentoc.py	(revision 5041682338ebf95a462e3296b92e84d6402011ff)
+++ doc/gentoc.py	(revision a5b8477c53695e7df03ecec7346f0c1bd770d3f7)
@@ -9,4 +9,10 @@
 from sasmodels.core import load_model_info
 
+try:
+    from typing import Optional, BinaryIO, List, Dict
+except ImportError:
+    pass
+else:
+    from sasmodels.modelinfo import ModelInfo
 
 TEMPLATE="""\
@@ -27,4 +33,5 @@
 
 def _make_category(category_name, label, title, parent=None):
+    # type: (str, str, str, Optional[BinaryIO]) -> BinaryIO
     file = open(joinpath(MODEL_TOC_PATH, category_name+".rst"), "w")
     file.write(TEMPLATE%{'label':label, 'title':title, 'bar':'*'*len(title)})
@@ -34,10 +41,13 @@
 
 def _add_subcategory(category_name, parent):
+    # type: (str, BinaryIO) -> None
     parent.write("    %s.rst\n"%category_name)
 
 def _add_model(file, model_name):
+    # type: (IO[str], str) -> None
     file.write("    ../../model/%s.rst\n"%model_name)
 
 def _maybe_make_category(category, models, cat_files, model_toc):
+    # type: (str, List[str], Dict[str, BinaryIO], BinaryIO) -> None
     if category not in cat_files:
         print("Unexpected category %s containing"%category, models, file=sys.stderr)
@@ -46,9 +56,10 @@
 
 def generate_toc(model_files):
+    # type: (List[str]) -> None
     if not model_files:
         print("gentoc needs a list of model files", file=sys.stderr)
 
     # find all categories
-    category = {}
+    category = {} # type: Dict[str, List[str]]
     for item in model_files:
         # assume model is in sasmodels/models/name.py, and ignore the full path
@@ -56,8 +67,8 @@
         if model_name.startswith('_'): continue
         model_info = load_model_info(model_name)
-        if model_info['category'] is None:
+        if model_info.category is None:
             print("Missing category for", item, file=sys.stderr)
         else:
-            category.setdefault(model_info['category'],[]).append(model_name)
+            category.setdefault(model_info.category,[]).append(model_name)
 
     # Check category names
