Changeset 5b0335b in sasmodels

Timestamp:

Apr 4, 2016 9:36:11 AM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

68e7f9d

Parents:

885ee37 (diff), 204fd9b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge remote-tracking branch 'origin/master' into polydisp

Conflicts:

sasmodels/core.py
sasmodels/models/spherical_sld.py

Files:

• sasmodels/core.py (modified) (5 diffs)
• sasmodels/models/core_multi_shell.py (modified) (5 diffs)
• sasmodels/models/raspberry.py (modified) (2 diffs)
• sasmodels/models/spherical_sld.py (modified) (2 diffs)
• doc/developer/calculator.rst (added)
• doc/developer/index.rst (modified) (1 diff)
• example/polynomial.py (added)
• sasmodels/bumps_model.py (modified) (1 diff)
• sasmodels/compare.py (modified) (5 diffs)
• sasmodels/convert.py (modified) (1 diff)
• sasmodels/direct_model.py (modified) (2 diffs)
• sasmodels/generate.py (modified) (17 diffs)
• sasmodels/kernel_header.c (added)
• sasmodels/kernel_iq.c (added)
• sasmodels/kernelcl.py (modified) (9 diffs)
• sasmodels/kerneldll.py (modified) (9 diffs)
• sasmodels/kernelpy.py (modified) (2 diffs)
• sasmodels/mixture.py (modified) (4 diffs)
• sasmodels/model_test.py (modified) (2 diffs)
• sasmodels/modelinfo.py (added)
• sasmodels/models/cylinder.c (modified) (3 diffs)
• sasmodels/models/gel_fit.c (modified) (4 diffs)
• sasmodels/models/lib/Si.c (modified) (1 diff)
• sasmodels/models/lib/polevl.c (modified) (2 diffs)
• sasmodels/models/lib/sas_J0.c (modified) (2 diffs)
• sasmodels/models/lib/sas_J1.c (modified) (2 diffs)
• sasmodels/models/lib/sas_JN.c (modified) (1 diff)
• sasmodels/models/lib/sph_j1c.c (modified) (2 diffs)
• sasmodels/models/lib/sphere_form.c (modified) (1 diff)
• sasmodels/models/lib/wrc_cyl.c (modified) (2 diffs)
• sasmodels/models/onion.c (modified) (6 diffs)
• sasmodels/models/onion.py (modified) (3 diffs)
• sasmodels/models/rpa.py (modified) (1 diff)
• sasmodels/product.py (modified) (3 diffs)
• sasmodels/resolution.py (modified) (5 diffs)
• sasmodels/sasview_model.py (modified) (10 diffs)
• sasmodels/src (added)

sasmodels/core.py

@@ -2 +2 @@
 Core model handling routines.
 """
+__all__ = [
+    "list_models", "load_model_info", "precompile_dll",
+    "build_model", "make_kernel", "call_kernel", "call_ER_VR",
+    ]
 
 from os.path import basename, dirname, join as joinpath, splitext
@@ -39 +43 @@
         return module
 
-
-
-__all__ = [
-    "list_models", "load_model_info", "precompile_dll",
-    "build_model", "call_kernel", "call_ER_VR",
-]
+try:
+    np.meshgrid([])
+    meshgrid = np.meshgrid
+except ValueError:
+    # CRUFT: np.meshgrid requires multiple vectors
+    def meshgrid(*args):
+        if len(args) > 1:
+            return np.meshgrid(*args)
+        else:
+            return [np.asarray(v) for v in args]
 
 def list_models():
@@ -90 +98 @@
         return product.make_product_info(P_info, Q_info)
 
-    return make_model_by_name(model_name)
-
-
-def make_model_by_name(model_name):
+    kernel_module = load_kernel_module(model_name)
+    return generate.make_model_info(kernel_module)
+
+
+def load_kernel_module(model_name):
     if model_name.endswith('.py'):
         path = model_name
@@ -106 +115 @@
         kernel_module = getattr(models, model_name, None)
     #import sys; print "\n".join(sys.path)
-    return generate.make_model_info(kernel_module)
+    __import__('sasmodels.models.'+model_name)
+    kernel_module = getattr(models, model_name, None)
+    return kernel_module
 
 
@@ -210 +221 @@
     """
     value, weight = zip(*pars)
-    value = [v.flatten() for v in np.meshgrid(*value)]
-    weight = np.vstack([v.flatten() for v in np.meshgrid(*weight)])
+    value = [v.flatten() for v in meshgrid(*value)]
+    weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
     weight = np.prod(weight, axis=0)
     return value, weight

sasmodels/models/core_multi_shell.py

@@ -89 +89 @@
 parameters = [["volfraction", "", 0.05, [0,1],"",
                "volume fraction of spheres"],
-              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+              ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
                "Core scattering length density"],
-              ["radius", "Ang", 200., [0, inf], "",
+              ["radius", "Ang", 200., [0, inf], "volume",
                "Radius of the core"],
               ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
@@ -97 +97 @@
               ["n", "", 1, [0, 10], "volume",
                "number of shells"],
-              ["sld_shell[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
+              ["sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
                "scattering length density of shell k"],
-              ["thick_shell[n]", "Ang", 40., [0, inf], "volume",
+              ["thickness[n]", "Ang", 40., [0, inf], "volume",
                "Thickness of shell k"],
               ]
@@ -112 +112 @@
 
 
-def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A):
+def profile(volfraction, sld_core, radius, sld_solvent, n, sld, thickness):
     """
     SLD profile
@@ -153 +153 @@
 # cut an the onion.  Seems like we should be consistant?
 
-    total_radius = 1.25*(sum(thickness[:n]) + core_radius + 1)
-    dr = total_radius/400  # 400 points for a smooth plot
+    total_radius = 1.25*(sum(thickness[:n]) + radius + 1)
 
     r = []
@@ -170 +169 @@
         r0 = r[-1]
         r.append(r0)
-        beta.append(sld_shell[k])
+        beta.append(sld[k])
         r.append(r0 + thickness[k])
-        beta.append(sld_shell[k])
+        beta.append(sld[k])
    # add in the solvent
     r.append(r[-1])
-    beta.append(solvent_sld)
+    beta.append(sld_solvent)
     r.append(total_radius)
-    beta.append(solvent_sld)
+    beta.append(sld_solvent)
 
     return np.asarray(r), np.asarray(beta)
 
-def ER(radius, n, thick_shell):
-    return np.sum(thick_shell[:n], axis=0) + core_radius
+def ER(radius, n, thickness):
+    n = n[0]  # n cannot be polydisperse
+    return np.sum(thickness[:n], axis=0) + radius
 
 def VR(radius, n, thick_shell):
-    return 1.0
-
-parameters = [["volfraction", "", 0.05, [0,1],"",
-               "volume fraction of spheres"],
-              ["sld", "1e-6/Ang^2", 1.0, [-inf, inf], "",
-               "Core scattering length density"],
-              ["radius", "Ang", 200., [0, inf], "",
-               "Radius of the core"],
-              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
-               "Solvent scattering length density"],
-              ["n", "", 1, [0, 10], "volume",
-               "number of shells"],
-              ["sld_shell[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
-               "scattering length density of shell k"],
-              ["thick_shell[n]", "Ang", 40., [0, inf], "volume",
-               "Thickness of shell k"],
-               ]
+    return 1.0, 1.0
 
 demo = dict(volfraction = 1.0,

sasmodels/models/raspberry.py

@@ -129 +129 @@
             Ref: J. coll. inter. sci. (2010) vol. 343 (1) pp. 36-41."""
 category = "shape:sphere"
+#single = False
 
 #             [ "name", "units", default, [lower, upper], "type", "description"],
@@ -163 +164 @@
 # names and the target sasview model name.
 
+# TODO: update tests so the parameters correspond to SasView parameters
+# The model was re-parameterized so the results have changed.
 # NOTE: test results taken from values returned by SasView 3.1.2, with
 # 0.001 added for a non-zero default background.
-tests = [[{}, 0.0412755102041, 0.286669115234],
-         [{}, 0.5, 0.00103818393658],
-        ]
+#tests = [[{}, 0.0412755102041, 0.286669115234],
+#         [{}, 0.5, 0.00103818393658],
+#        ]

sasmodels/models/spherical_sld.py

@@ -168 +168 @@
 category = "sphere-based"
 
+INTERFACE_CASES = ["Erf", "RPower", "LPower", "RExp", "LExp"]
+
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
 parameters = [["n_shells",                "",               1,      [0, 9],      "", "number of shells"],
-              ["thick_inter[n_shells]",   "Ang",            50,     [-inf, inf], "", "intern layer thickness"],
-              ["func_inter[n_shells]",    "",               0,      [0, 4],      "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4"],
+              ["radius_core",             "Ang",            50.0,   [0, inf],    "", "intern layer thickness"],
               ["sld_core",                "1e-6/Ang^2",     2.07,   [-inf, inf], "", "sld function flat"],
               ["sld_solvent",             "1e-6/Ang^2",     1.0,    [-inf, inf], "", "sld function solvent"],
               ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf], "", "sld function flat"],
+              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],    "", "flat layer_thickness"],
+              ["func_inter[n_shells]",    "",               0,      INTERFACE_CASES,  "", "shape of the interface between shells"],
+              ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf], "", "steepness parameter"],
               ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],    "", "intern layer thickness"],
-              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],    "", "flat layer_thickness"],
-              ["inter_nu[n_shells]",      "",               2.5,    [-inf, inf], "", "steepness parameter"],
-              ["npts_inter",              "",               35,     [0, 35],     "", "number of points in each sublayer Must be odd number"],
-              ["core_rad",                "Ang",            50.0,   [0, inf],    "", "intern layer thickness"],
+              ["n_points_inter",          "",               35,     [0, 35],     "", "number of points in each sublayer Must be odd number"],
               ]
 # pylint: enable=bad-whitespace, line-too-long
@@ -193 +194 @@
     solvent_sld=1.0,
     background=0.0,
-    npts_inter=35.0,
-    func_inter_0=0,
-    nu_inter_0=2.5,
-    rad_core_0=50.0,
-    core0_sld=2.07,
-    thick_inter_0=50.0,
-    func_inter_1=0,
-    nu_inter_1=2.5,
-    thick_inter_1=50.0,
-    flat1_sld=4.0,
-    thick_flat_1=100.0,
-    func_inter_2=0,
-    nu_inter_2=2.5,
-    thick_inter_2=50.0,
-    flat2_sld=3.5,
-    thick_flat_2=100.0,
-    func_inter_3=0,
-    nu_inter_3=2.5,
-    thick_inter_3=50.0,
-    flat3_sld=4.0,
-    thick_flat_3=100.0,
-    func_inter_4=0,
-    nu_inter_4=2.5,
-    thick_inter_4=50.0,
-    flat4_sld=3.5,
-    thick_flat_4=100.0,
-    func_inter_5=0,
-    nu_inter_5=2.5,
-    thick_inter_5=50.0,
-    flat5_sld=4.0,
-    thick_flat_5=100.0,
-    func_inter_6=0,
-    nu_inter_6=2.5,
-    thick_inter_6=50.0,
-    flat6_sld=3.5,
-    thick_flat_6=100.0,
-    func_inter_7=0,
-    nu_inter_7=2.5,
-    thick_inter_7=50.0,
-    flat7_sld=4.0,
-    thick_flat_7=100.0,
-    func_inter_8=0,
-    nu_inter_8=2.5,
-    thick_inter_8=50.0,
-    flat8_sld=3.5,
-    thick_flat_8=100.0,
-    func_inter_9=0,
-    nu_inter_9=2.5,
-    thick_inter_9=50.0,
-    flat9_sld=4.0,
-    thick_flat_9=100.0,
-    func_inter_10=0,
-    nu_inter_10=2.5,
-    thick_inter_10=50.0,
-    flat10_sld=3.5,
-    thick_flat_10=100.0
+    n_points_inter=35.0,
     )
 
 oldname = "SphereSLDModel"
 oldpars = dict(
-    n_shells="n_shells",
-    scale="scale",
-    background='background',
-    npts_inter='npts_inter',
-    solvent_sld='sld_solv',
-
-    func_inter_0='func_inter0',
-    nu_inter_0='nu_inter0',
-    rad_core_0='rad_core0',
-    core0_sld='sld_core0',
-    thick_inter_0='thick_inter0',
-    func_inter_1='func_inter1',
-    nu_inter_1='nu_inter1',
-    thick_inter_1='thick_inter1',
-    flat1_sld='sld_flat1',
-    thick_flat_1='thick_flat1',
-    func_inter_2='func_inter2',
-    nu_inter_2='nu_inter2',
-    thick_inter_2='thick_inter2',
-    flat2_sld='sld_flat2',
-    thick_flat_2='thick_flat2',
-    func_inter_3='func_inter3',
-    nu_inter_3='nu_inter3',
-    thick_inter_3='thick_inter3',
-    flat3_sld='sld_flat3',
-    thick_flat_3='thick_flat3',
-    func_inter_4='func_inter4',
-    nu_inter_4='nu_inter4',
-    thick_inter_4='thick_inter4',
-    flat4_sld='sld_flat4',
-    thick_flat_4='thick_flat4',
-    func_inter_5='func_inter5',
-    nu_inter_5='nu_inter5',
-    thick_inter_5='thick_inter5',
-    flat5_sld='sld_flat5',
-    thick_flat_5='thick_flat5',
-    func_inter_6='func_inter6',
-    nu_inter_6='nu_inter6',
-    thick_inter_6='thick_inter6',
-    flat6_sld='sld_flat6',
-    thick_flat_6='thick_flat6',
-    func_inter_7='func_inter7',
-    nu_inter_7='nu_inter7',
-    thick_inter_7='thick_inter7',
-    flat7_sld='sld_flat7',
-    thick_flat_7='thick_flat7',
-    func_inter_8='func_inter8',
-    nu_inter_8='nu_inter8',
-    thick_inter_8='thick_inter8',
-    flat8_sld='sld_flat8',
-    thick_flat_8='thick_flat8',
-    func_inter_9='func_inter9',
-    nu_inter_9='nu_inter9',
-    thick_inter_9='thick_inter9',
-    flat9_sld='sld_flat9',
-    thick_flat_9='thick_flat9',
-    func_inter_10='func_inter10',
-    nu_inter_10='nu_inter10',
-    thick_inter_10='thick_inter10',
-    flat10_sld='sld_flat10',
-    thick_flat_10='thick_flat10')
+    #scale="scale",
+    #background='background',
+    #n_shells="n_shells",
+    radius_core='rad_core',
+    #sld_core='sld_core',
+    #sld_flat='sld_flat',
+    #thick_flat='thick_flat',
+    #func_inter='func_inter',
+    #thick_inter='thick_inter',
+    #nu_inter='nu_inter',
+    #points_inter='npts_inter',
+    sld_solvent='sld_solv',
+    )
 
 #TODO: Not working yet
 tests = [
     # Accuracy tests based on content in test/utest_extra_models.py
-    [{'npts_iter':35,
-        'sld_solv':1,
-        'func_inter_1':0.0,
-        'nu_inter':2.5,
-        'thick_inter_1':50,
-        'sld_flat_1':4,
-        'thick_flat_1':100,
-        'func_inter_0':0.0,
-        'nu_inter_0':2.5,
-        'rad_core_0':50.0,
-        'sld_core_0':2.07,
-        'thick_inter_0':50,
-        'background': 0.0,
-    }, 0.001, 1000],
+    [{'n_points_inter':35,
+      'sld_solv':1,
+      'radius_core':50.0,
+      'sld_core':2.07,
+      'func_inter2':0.0,
+      'thick_inter2':50,
+      'nu_inter2':2.5,
+      'sld_flat2':4,
+      'thick_flat2':100,
+      'func_inter1':0.0,
+      'thick_inter1':50,
+      'nu_inter1':2.5,
+      'background': 0.0,
+     }, 0.001, 0.001],
 ]

doc/developer/index.rst

@@ -3 +3 @@
 ***************************
 
+.. toctree::
+   :numbered: 4
+   :maxdepth: 4
 
+   calculator.rst

sasmodels/bumps_model.py

@@ -81 +81 @@
     from bumps.names import Parameter
 
-    pars = {}
+    pars = {}     # floating point parameters
+    pd_types = {} # distribution names
     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 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)
+            pd_types[p.name+'_pd_type'] = kwargs.pop(name, 'gaussian')
 
     if kwargs:  # args not corresponding to parameters

sasmodels/compare.py

@@ -306 +306 @@
     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)
@@ -454 +448 @@
     """
     # 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
@@ -661 +656 @@
     """
     # 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
@@ -774 +773 @@
 
     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')
@@ -880 +879 @@
         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 p in model_info['parameters'].type['1d']:
+                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/convert.py

@@ -73 +73 @@
     new model definition end with sld.
     """
+    print "pars",pars
     return dict((p, (v*1e-6 if p.startswith('sld') or p.endswith('sld')
                      else v*1e15 if 'ndensity' in p

sasmodels/direct_model.py

@@ -68 +68 @@
             self.data_type = 'Iq'
 
-        partype = model.info['partype']
-
         if self.data_type == 'sesans':
             
@@ -83 +81 @@
             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)

sasmodels/generate.py

@@ -21 +21 @@
 
     *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
@@ -65 +69 @@
 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:
 
@@ -212 +181 @@
 from __future__ import print_function
 
-# TODO: identify model files which have changed since loading and reload them.
-
-import sys
+#TODO: determine which functions are useful outside of generate
+#__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
+
 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
 
-PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description']
-Parameter = namedtuple('Parameter', PARAMETER_FIELDS)
-
-#TODO: determine which functions are useful outside of generate
-#__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
-
-C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
+from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo
+
+# TODO: identify model files which have changed since loading and reload them.
+
+TEMPLATE_ROOT = dirname(__file__)
 
 F16 = np.dtype('float16')
@@ -239 +205 @@
 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
@@ -338 +298 @@
     raise ValueError("%r not found in %s" % (filename, search_path))
 
+
 def model_sources(model_info):
     """
@@ -346 +307 @@
     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):
@@ -369 +327 @@
     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 +333 @@
     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 +376 @@
 
 
-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):
     """
@@ -460 +462 @@
     # 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_coord(self): return self.details[-2]
+    @num_coord.setter
+    def num_coord(self, v): self.details[-2] = 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_active(self): return self.details[-4]
+    @num_active.setter
+    def num_active(self, v): self.details[-4] = v
+
+    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_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 callable(model_info['Iq']) and model_info['Iqxy'] is None:
+        partable = model_info['parameters']
+        if partable.has_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):
@@ -678 +689 @@
       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
@@ -699 +706 @@
       *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]
@@ -712 +719 @@
         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,
@@ -720 +727 @@
         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),
@@ -727 +735 @@
         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
 
@@ -782 +795 @@
 
 
-
 def demo_time():
     """
@@ -798 +810 @@
     Program which prints the source produced by the model.
     """
+    import sys
+    from sasmodels.core import make_model_by_name
     if len(sys.argv) <= 1:
         print("usage: python -m sasmodels.generate modelname")
     else:
         name = sys.argv[1]
-        import sasmodels.models
-        __import__('sasmodels.models.' + name)
-        model = getattr(sasmodels.models, name)
-        model_info = make_model_info(model)
+        model_info = make_model_by_name(name)
         source = make_source(model_info)
         print(source)

sasmodels/kernelcl.py

@@ -48 +48 @@
 harmless, albeit annoying.
 """
+from __future__ import print_function
 import os
 import warnings
@@ -73 +74 @@
 # 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
+"""
 
 
@@ -142 +158 @@
         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
@@ -220 +243 @@
         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]
@@ -314 +338 @@
         self.program = None
 
-    def __call__(self, q_vectors):
+    def make_kernel(self, q_vectors):
         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)
@@ -365 +389 @@
         # 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):
@@ -378 +412 @@
         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):
@@ -406 +440 @@
     """
     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
@@ -418 +455 @@
         env = environment()
         self.queue = env.get_queue(dtype)
-        self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                                 2 * 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

@@ -50 +50 @@
 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
@@ -81 +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 = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm"
@@ -90 +89 @@
 
 
-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"):
@@ -133 +132 @@
         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}
@@ -171 +162 @@
     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):
     """
@@ -197 +184 @@
 
     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:
@@ -212 +194 @@
               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):
@@ -229 +209 @@
         self.dll = None
 
-    def __call__(self, q_vectors):
+    def make_kernel(self, q_vectors):
         q_input = PyInput(q_vectors, self.dtype)
         if self.dll is None: self._load_dll()
@@ -272 +252 @@
     """
     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

@@ -53 +53 @@
         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):
@@ -60 +65 @@
         Free resources associated with the model inputs.
         """
-        self.q_vectors = []
+        self.q = None
 
 class PyKernel(object):

sasmodels/mixture.py

@@ -14 +14 @@
 import numpy as np
 
-from .generate import process_parameters, COMMON_PARAMETERS, Parameter
+from .modelinfo import Parameter, ParameterTable
 
 SCALE=0
@@ -34 +34 @@
 
     # 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 = {}
@@ -58 +58 @@
     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
@@ -67 +67 @@
     # Remember the component info blocks so we can build the model
     model_info['composition'] = ('mixture', parts)
-    process_parameters(model_info)
-    return model_info
 
 

sasmodels/model_test.py

@@ -51 +51 @@
 
 from .core import list_models, load_model_info, build_model, HAVE_OPENCL
-from .core import make_kernel, call_kernel, call_ER, call_VR
+from .core import call_kernel, call_ER, call_VR
 from .exception import annotate_exception
 
@@ -187 +187 @@
                 Qx, Qy = zip(*x)
                 q_vectors = [np.array(Qx), np.array(Qy)]
-                kernel = make_kernel(model, q_vectors)
+                kernel = model.make_kernel(q_vectors)
                 actual = call_kernel(kernel, pars)
             else:
                 q_vectors = [np.array(x)]
-                kernel = make_kernel(model, q_vectors)
+                kernel = model.make_kernel(q_vectors)
                 actual = call_kernel(kernel, pars)
 

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/models/lib/Si.c

@@ +1 @@
+// integral of sin(x)/x Taylor series approximated to w/i 0.1%
 double Si(double x);
 
-// integral of sin(x)/x Taylor series approximated to w/i 0.1%
 double Si(double x)
 {

sasmodels/models/lib/polevl.c

@@ -51 +51 @@
 */
 
-double polevl( double x, constant double *coef, int N ) {
+double polevl( double x, constant double *coef, int N );
+double p1evl( double x, constant double *coef, int N );
+
+
+double polevl( double x, constant double *coef, int N )
+{
 
     int i = 0;
@@ -70 +75 @@
  */
 
-double p1evl( double x, constant double *coef, int N ) {
+double p1evl( double x, constant double *coef, int N )
+{
     int i=0;
     double ans = x+coef[i];

sasmodels/models/lib/sas_J0.c

@@ -49 +49 @@
 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
 */
+
+double sas_J0(double x);
 
 /* Note: all coefficients satisfy the relative error criterion
@@ -177 +179 @@
  };
 
-double sas_J0(double x) {
+double sas_J0(double x)
+{
 
 //Cephes single precission

sasmodels/models/lib/sas_J1.c

@@ -39 +39 @@
 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
 */
+double sas_J1(double x);
+double sas_J1c(double x);
 
 constant double RPJ1[8] = {
@@ -135 +137 @@
     };
 
-double sas_J1(double x) {
+double sas_J1(double x)
+{
 
 //Cephes double pression function

sasmodels/models/lib/sas_JN.c

@@ -47 +47 @@
 Copyright 1984, 1987, 2000 by Stephen L. Moshier
 */
-
 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

@@ -7 +7 @@
 * 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
@@ -43 +44 @@
 #endif
 
-double sph_j1c(double q);
 double sph_j1c(double q)
 {

sasmodels/models/lib/sphere_form.c

@@ -1 @@
-inline double
-sphere_volume(double radius)
+double sphere_volume(double radius);
+double sphere_form(double q, double radius, double sld, double solvent_sld);
+
+double sphere_volume(double radius)
 {
     return M_4PI_3*cube(radius);
 }
 
-inline double
-sphere_form(double q, double radius, double sld, double solvent_sld)
+double sphere_form(double q, double radius, double sld, double solvent_sld)
 {
     const double fq = sphere_volume(radius) * sph_j1c(q*radius);

sasmodels/models/lib/wrc_cyl.c

@@ -2 +2 @@
     Functions for WRC implementation of flexible cylinders
 */
+double Sk_WR(double q, double L, double b);
+
+
 static double
 AlphaSquare(double x)
@@ -363 +366 @@
 }
 
-double Sk_WR(double q, double L, double b);
 double Sk_WR(double q, double L, double b)
 {

sasmodels/models/onion.c

@@ -4 +4 @@
     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;
@@ -19 +19 @@
     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)
@@ -32 +32 @@
 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);
@@ -52 +52 @@
 {
   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;
 }
@@ -64 +64 @@
     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;
@@ -92 +92 @@
     }
   }
-  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

@@ -293 +293 @@
 
 #             ["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",
@@ -311 +311 @@
               ]
 
-#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
@@ -374 +371 @@
 
 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.py

@@ -86 +86 @@
 #   ["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/product.py

@@ -14 +14 @@
 
 from .core import call_ER_VR
-from .generate import process_parameters
 
 SCALE=0
@@ -66 +65 @@
     model_info['oldpars'] = oldpars
     model_info['composition'] = ('product', [p_info, s_info])
-    process_parameters(model_info)
     return model_info
 
@@ -101 +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

@@ -477 +477 @@
     """
     from sasmodels import core
-    kernel = core.make_kernel(form, [q])
+    kernel = form.make_kernel([q])
     theory = core.call_kernel(kernel, pars)
     kernel.release()
@@ -502 +502 @@
     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))))
@@ -556 +556 @@
     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))))
+                         (", ".join(sorted(extra)),
+                          ", ".join(sorted(pars.keys()))))
 
     _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
@@ -693 +694 @@
     def _eval_sphere(self, pars, resolution):
         from sasmodels import core
-        kernel = core.make_kernel(self.model, [resolution.q_calc])
+        kernel = self.model.make_kernel([resolution.q_calc])
         theory = core.call_kernel(kernel, pars)
         result = resolution.apply(theory)
@@ -1062 +1063 @@
     model = core.build_model(model_info)
 
-    kernel = core.make_kernel(model, [resolution.q_calc])
+    kernel = model.make_kernel([resolution.q_calc])
     theory = core.call_kernel(kernel, pars)
     Iq = resolution.apply(theory)

sasmodels/sasview_model.py

@@ -21 +21 @@
 
 from . import core
-from . import generate
+from . import weights
 
 def standard_models():
@@ -32 +32 @@
 
     Returns a class that can be used directly as a sasview model.t
-
-    Defaults to using the new name for a model.  Setting
-    *namestyle='oldname'* will produce a class with a name
-    compatible with SasView.
     """
     model_info = core.load_model_info(model_name)
@@ -56 +52 @@
     """
     def __init__(self):
-        self._kernel = None
+        self._model = None
         model_info = self._model_info
+        parameters = model_info['parameters']
 
         self.name = model_info['name']
-        self.oldname = model_info['oldname']
         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
@@ -71 +75 @@
         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 = []
 
@@ -109 +116 @@
     def __get_state__(self):
         state = self.__dict__.copy()
-        model_id = self._model_info['id']
         state.pop('_kernel')
         # May need to reload model info on set state since it has pointers
@@ -118 +124 @@
     def __set_state__(self, state):
         self.__dict__ = state
-        self._kernel = None
+        self._model = None
 
     def __str__(self):
@@ -207 +213 @@
     def getDispParamList(self):
         """
-        Return a list of all available parameters for the model
+        Return a list of polydispersity parameters for the model
         """
         # 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
@@ -285 +292 @@
             # 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:
@@ -308 +314 @@
         to the card for each evaluation.
         """
-        if self._kernel is None:
-            self._kernel = core.build_model(self._model_info)
+        if self._model is None:
+            self._model = core.build_model(self._model_info, platform='dll')
         q_vectors = [np.asarray(q) for q in args]
-        fn = self._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)
+        return kernel(details, weights, values, cutoff=self.cutoff)
+        kernel.q_input.release()
+        kernel.release()
         return result
 
@@ -389 +396 @@
     def _get_weights(self, par):
         """
-            Return dispersion weights
-            :param par parameter name
-        """
-        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)
-
+        Return dispersion weights for parameter
+        """
+        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():
+    Cylinder = make_class('cylinder')
+    cylinder = Cylinder()
+    return cylinder.evalDistribution([0.1,0.1])
+
+if __name__ == "__main__":
+    print("cylinder(0.1,0.1)=%g"%test_model())