Changeset 1e2a1ba in sasmodels

Timestamp:

Apr 5, 2016 12:06:36 PM (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:

d2bb604

Parents:

ea05c87 (diff), 4d76711 (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/custom/init.py
sasmodels/direct_model.py
sasmodels/generate.py
sasmodels/kernelpy.py
sasmodels/sasview_model.py

Files:

• sasmodels/core.py (modified) (2 diffs)
• sasmodels/custom/__init__.py (modified) (1 diff)
• sasmodels/exception.py (modified) (4 diffs)
• sasmodels/generate.py (modified) (2 diffs)
• sasmodels/kerneldll.py (modified) (2 diffs)
• sasmodels/kernelpy.py (modified) (1 diff)
• sasmodels/model_test.py (modified) (1 diff)
• sasmodels/sasview_model.py (modified) (4 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/direct_model.py (modified) (3 diffs)
• sasmodels/kernel_header.c (added)
• sasmodels/kernel_iq.c (added)
• sasmodels/kernelcl.py (modified) (9 diffs)
• sasmodels/mixture.py (modified) (4 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) (2 diffs)
• sasmodels/src (added)

sasmodels/core.py

@@ -38 +38 @@
             return [np.asarray(v) for v in args]
 
+# TODO: refactor composite model support
+# The current load_model_info/build_model does not reuse existing model
+# definitions when loading a composite model, instead reloading and
+# rebuilding the kernel for each component model in the expression.  This
+# is fine in a scripting environment where the model is built when the script
+# starts and is thrown away when the script ends, but may not be the best
+# solution in a long-lived application.  This affects the following functions:
+#
+#    load_model
+#    load_model_info
+#    build_model
+
 def list_models():
     """
@@ -60 +72 @@
     """
     return build_model(load_model_info(model_name), **kw)
+
 
 def load_model_info(model_name):

sasmodels/custom/__init__.py

@@ -32 +32 @@
     name = basename(splitext(path)[0])
     # Placing the model in the 'sasmodels.custom' name space.
-    from sasmodels import custom
     kernel_module = load_module_from_path('sasmodels.custom.'+name, path)
     return kernel_module

sasmodels/exception.py

@@ -2 +2 @@
 Utility to add annotations to python exceptions.
 """
+import sys
 
 # Platform cruft: WindowsError is only defined on Windows.
@@ -13 +14 @@
         pass
 
-def annotate_exception(exc, msg):
+def annotate_exception(msg, exc=None):
     """
     Add an annotation to the current exception, which can then be forwarded
-    to the caller using a bare "raise" statement to reraise the annotated
-    exception.
+    to the caller using a bare "raise" statement to raise the annotated
+    exception.  If the exception *exc* is provided, then that exception is the
+    one that is annotated, otherwise *sys.exc_info* is used.
 
     Example::
@@ -24 +26 @@
         >>> try:
         ...    print(D['hello'])
-        ... except Exception as exc:
-        ...    annotate_exception(exc, "while accessing 'D'")
+        ... except:
+        ...    annotate_exception("while accessing 'D'")
         ...    raise
         Traceback (most recent call last):
@@ -31 +33 @@
         KeyError: "hello while accessing 'D'"
     """
+    if not exc:
+        exc = sys.exc_info()[1]
+
     # Can't extend WindowsError exceptions; instead raise a new exception.
     # TODO: try to incorporate current stack trace in the raised exception

sasmodels/generate.py

@@ -186 +186 @@
 from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo
 from .custom import load_custom_kernel_module
-
-# TODO: identify model files which have changed since loading and reload them.
 
 TEMPLATE_ROOT = dirname(__file__)
@@ -320 +318 @@
     if dtype == F16:
         fbytes = 2
-        source = _convert_type(source, "float", "f")
+        source = _convert_type(source, "half", "f")
     elif dtype == F32:
         fbytes = 4

sasmodels/kerneldll.py

@@ -187 +187 @@
         try:
             self.dll = ct.CDLL(self.dllpath)
-        except Exception as exc:
-            annotate_exception(exc, "while loading "+self.dllpath)
+        except:
+            annotate_exception("while loading "+self.dllpath)
             raise
 
@@ -214 +214 @@
         kernel = self.Iqxy if q_input.is_2d else self.Iq
         return DllKernel(kernel, self.info, q_input)
-    
+
     def release(self):
         """

sasmodels/kernelpy.py

@@ -19 +19 @@
         self.info = model_info
 
-    def __call__(self, q_vectors):
+    def make_kernel(self, q_vectors):
         q_input = PyInput(q_vectors, dtype=F64)
         kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq']

sasmodels/model_test.py

@@ -166 +166 @@
                     pass
 
-            except Exception as exc:
-                annotate_exception(exc, self.test_name)
+            except:
+                annotate_exception(self.test_name)
                 raise
 

sasmodels/sasview_model.py

@@ -13 +13 @@
 using :func:`sasmodels.convert.convert`.
 """
+from __future__ import print_function
 
 import math
 from copy import deepcopy
 import collections
+import traceback
+import logging
 
 import numpy as np
 
 from . import core
+from . import custom
+from . import generate
 from . import weights
 
-def standard_models():
-    return [make_class(model_name) for model_name in core.list_models()]
-
-# TODO: rename to make_class_from_name and update sasview
-def make_class(model_name):
-    """
-    Load the sasview model defined in *kernel_module*.
-
-    Returns a class that can be used directly as a sasview model.t
-    """
-    model_info = core.load_model_info(model_name)
-    return make_class_from_info(model_info)
-
-def make_class_from_file(path):
-    model_info = core.load_model_info_from_path(path)
-    return make_class_from_info(model_info)
-
-def make_class_from_info(model_info):
+def load_standard_models():
+    """
+    Load and return the list of predefined models.
+
+    If there is an error loading a model, then a traceback is logged and the
+    model is not returned.
+    """
+    models = []
+    for name in core.list_models():
+        try:
+            models.append(_make_standard_model(name))
+        except:
+            logging.error(traceback.format_exc())
+    return models
+
+
+def load_custom_model(path):
+    """
+    Load a custom model given the model path.
+    """
+    kernel_module = custom.load_custom_kernel_module(path)
+    model_info = generate.make_model_info(kernel_module)
+    return _make_model_from_info(model_info)
+
+
+def _make_standard_model(name):
+    """
+    Load the sasview model defined by *name*.
+
+    *name* can be a standard model name or a path to a custom model.
+
+    Returns a class that can be used directly as a sasview model.
+    """
+    kernel_module = generate.load_kernel_module(name)
+    model_info = generate.make_model_info(kernel_module)
+    return _make_model_from_info(model_info)
+
+
+def _make_model_from_info(model_info):
+    """
+    Convert *model_info* into a SasView model wrapper.
+    """
     def __init__(self, multfactor=1):
         SasviewModel.__init__(self)
@@ -47 +76 @@
     return ConstructedModel
 
+
 class SasviewModel(object):
     """
     Sasview wrapper for opencl/ctypes model.
     """
+    _model_info = {}
     def __init__(self):
         self._model = None
@@ -116 +147 @@
     def __get_state__(self):
         state = self.__dict__.copy()
-        state.pop('_kernel')
+        state.pop('_model')
         # May need to reload model info on set state since it has pointers
         # to python implementations of Iq, etc.
@@ -408 +439 @@
 
 def test_model():
-    Cylinder = make_class('cylinder')
+    """
+    Test that a sasview model (cylinder) can be run.
+    """
+    Cylinder = _make_standard_model('cylinder')
     cylinder = Cylinder()
     return cylinder.evalDistribution([0.1,0.1])
 
+
+def test_model_list():
+    """
+    Make sure that all models build as sasview models.
+    """
+    from .exception import annotate_exception
+    for name in core.list_models():
+        try:
+            _make_standard_model(name)
+        except:
+            annotate_exception("when loading "+name)
+            raise
+
 if __name__ == "__main__":
     print("cylinder(0.1,0.1)=%g"%test_model())

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/direct_model.py

@@ -67 +67 @@
             self.data_type = 'Iq'
 
-        partype = model.info['partype']
-
         if self.data_type == 'sesans':
             
@@ -82 +80 @@
             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)
@@ -172 +170 @@
     def _calc_theory(self, pars, cutoff=0.0):
         if self._kernel is None:
-            self._kernel = self._model.make_kernel(self._kernel_inputs)  # pylint: disable=attribute-dedata_type
+            self._kernel = self._model.make_kernel(self._kernel_inputs)
             self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs)
                                  if self._kernel_mono_inputs else None)

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]
@@ -317 +341 @@
         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/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/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
@@ -58 +57 @@
     # Iq, Iqxy, form_volume, ER, VR and sesans
     model_info['composition'] = ('product', [p_info, s_info])
-    process_parameters(model_info)
     return model_info
 
@@ -93 +91 @@
         # 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

@@ -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)