Changeset 9150036 – SasView

doc/guide/plugin.rst

-                      raa8c6e0
+                      r9150036
 structure factor to account for interactions between particles.  See
 `Form_Factors`_ for more details.
+**model_info = ...** lets you define a model directly, for example, by
+loading and modifying existing models.  This is done implicitly by
+:func:`sasmodels.core.load_model_info`, which can create a mixture model
+from a pair of existing models.  For example::
+    from sasmodels.core import load_model_info
+    model_info = load_model_info('sphere+cylinder')
+See :class:`sasmodels.modelinfo.ModelInfo` for details about the model
+attributes that are defined.
 Model Parameters
 …
              - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)
         For small arguments ,
+        For small arguments,
         .. math::

example/multiscatfit.py

-                      r49d1f8b8
+                      r2c4a190
     # Show the model without fitting
     PYTHONPATH=..:../explore:../../bumps:../../sasview/src python multiscatfit.py
+    PYTHONPATH=..:../../bumps:../../sasview/src python multiscatfit.py
     # Run the fit
     PYTHONPATH=..:../explore:../../bumps:../../sasview/src ../../bumps/run.py \
+    PYTHONPATH=..:../../bumps:../../sasview/src ../../bumps/run.py \
     multiscatfit.py --store=/tmp/t1
 …
+    )
+# Tie the model to the data
+M = Experiment(data=data, model=model)
+# Stack mulitple scattering on top of the existing resolution function.
+M.resolution = MultipleScattering(resolution=M.resolution, probability=0.)
 # SET THE FITTING PARAMETERS
 model.radius_polar.range(15, 3000)
 …
 model.scale.range(0, 0.1)
+# Mulitple scattering probability parameter
+# HACK: the probability is stuffed in as an extra parameter to the experiment.
+probability = Parameter(name="probability", value=0.0)
+probability.range(0.0, 0.9)
+# The multiple scattering probability parameter is in the resolution function
+# instead of the scattering function, so access it through M.resolution
+M.scattering_probability.range(0.0, 0.9)
+M = Experiment(data=data, model=model, extra_pars={'probability': probability})
+# Stack mulitple scattering on top of the existing resolution function.
+# Because resolution functions in sasview don't have fitting parameters,
+# we instead allow the multiple scattering calculator to take a function
+# instead of a probability.  This function returns the current value of
+# the parameter. ** THIS IS TEMPORARY ** when multiple scattering is
+# properly integrated into sasmodels and sasview, its fittable parameter
+# will be treated like the model parameters.
+M.resolution = MultipleScattering(resolution=M.resolution,
+                                  probability=lambda: probability.value,
+                                  )
+M._kernel_inputs = M.resolution.q_calc
+# Let bumps know that we are fitting this experiment
 problem = FitProblem(M)

sasmodels/init.py

ra1ec908	r37f38ff
14	14	defining new models.
15	15	"""
16		__version__ = "0.98"
	16	__version__ = "0.99"
17	17
18	18	def data_files():

sasmodels/bumps_model.py

-                      r49d1f8b8
+                      r2c4a190
     # when bumps is not on the path.
     from bumps.names import Parameter # type: ignore
+    from bumps.parameter import Reference # type: ignore
 except ImportError:
     pass
 …
     def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None):
         # type: (Data, Model, float) -> None
+        # Allow resolution function to define fittable parameters.  We do this
+        # by creating reference parameters within the resolution object rather
+        # than modifying the object itself to use bumps parameters.  We need
+        # to reset the parameters each time the object has changed.  These
+        # additional parameters need to be returned from the fitting engine.
+        # To make them available to the user, they are added as top-level
+        # attributes to the experiment object.  The only change to the
+        # resolution function is that it needs an optional 'fittable' attribute
+        # which maps the internal name to the user visible name for the
+        # for the parameter.
+        self._resolution = None
+        self._resolution_pars = {}
         # remember inputs so we can inspect from outside
         self.name = data.filename if name is None else name
 …
         self._interpret_data(data, model.sasmodel)
         self._cache = {}
+        # CRUFT: no longer need extra parameters
+        # Multiple scattering probability is now retrieved directly from the
+        # multiple scattering resolution function.
         self.extra_pars = extra_pars
 …
         return len(self.Iq)
+    @property
+    def resolution(self):
+        return self._resolution
+    @resolution.setter
+    def resolution(self, value):
+        self._resolution = value
+        # Remove old resolution fitting parameters from experiment
+        for name in self._resolution_pars:
+            delattr(self, name)
+        # Create new resolution fitting parameters
+        res_pars = getattr(self._resolution, 'fittable', {})
+        self._resolution_pars = {
+            name: Reference(self._resolution, refname, name=name)
+            for refname, name in res_pars.items()
+        }
+        # Add new resolution fitting parameters as experiment attributes
+        for name, ref in self._resolution_pars.items():
+            setattr(self, name, ref)
     def parameters(self):
         # type: () -> Dict[str, Parameter]
 …
         """
         pars = self.model.parameters()
         if self.extra_pars:
+        if self.extra_pars is not None:
             pars.update(self.extra_pars)
+        pars.update(self._resolution_pars)
         return pars

sasmodels/direct_model.py

-                      rc1799d3
+                      r9150036
             else:
                 Iq, dIq = None, None
-            #self._theory = np.zeros_like(q)
-            q_vectors = [res.q_calc]
         elif self.data_type == 'Iqxy':
             #if not model.info.parameters.has_2d:
 …
             res = resolution2d.Pinhole2D(data=data, index=index,
                                          nsigma=3.0, accuracy=accuracy)
-            #self._theory = np.zeros_like(self.Iq)
-            q_vectors = res.q_calc
         elif self.data_type == 'Iq':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
 …
             else:
                 res = resolution.Perfect1D(data.x[index])
-            #self._theory = np.zeros_like(self.Iq)
-            q_vectors = [res.q_calc]
         elif self.data_type == 'Iq-oriented':
             index = (data.x >= data.qmin) & (data.x <= data.qmax)
 …
                                       qx_width=data.dxw[index],
                                       qy_width=data.dxl[index])
-            q_vectors = res.q_calc
         else:
             raise ValueError("Unknown data type") # never gets here
 …
         # Remember function inputs so we can delay loading the function and
         # so we can save/restore state
-        self._kernel_inputs = q_vectors
         self._kernel = None
         self.Iq, self.dIq, self.index = Iq, dIq, index
 …
         # type: (ParameterSet, float) -> np.ndarray
         if self._kernel is None:
+            self._kernel = self._model.make_kernel(self._kernel_inputs)
+            # TODO: change interfaces so that resolution returns kernel inputs
+            # Maybe have resolution always return a tuple, or maybe have
+            # make_kernel accept either an ndarray or a pair of ndarrays.
+            kernel_inputs = self.resolution.q_calc
+            if isinstance(kernel_inputs, np.ndarray):
+                kernel_inputs = (kernel_inputs,)
+            self._kernel = self._model.make_kernel(kernel_inputs)
         # Need to pull background out of resolution for multiple scattering

sasmodels/multiscat.py

-                      rb3703f5
+                      r2c4a190
     *probability* is related to the expected number of scattering
+    events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a
+    hack to allow probability to be a fitted parameter, the "value"
+    can be a function that takes no parameters and returns the current
+    value of the probability.  *coverage* determines how many scattering
+    steps to consider.  The default is 0.99, which sets $n$ such that
+    $1 \ldots n$ covers 99% of the Poisson probability mass function.
+    events in the sample $\lambda$ as $p = 1 - e^{-\lambda}$.
+    *coverage* determines how many scattering steps to consider.  The
+    default is 0.99, which sets $n$ such that $1 \ldots n$ covers 99%
+    of the Poisson probability mass function.
     *is2d* is True then 2D scattering is used, otherwise it accepts
 …
         self.qmin = qmin
         self.nq = nq
         self.probability = probability
+        self.probability = 0. if probability is None else probability
         self.coverage = coverage
         self.is2d = is2d
 …
         self.Iqxy = None # type: np.ndarray
+        # Label probability as a fittable parameter, and give its external name
+        # Note that the external name must be a valid python identifier, since
+        # is will be set as an experiment attribute.
+        self.fittable = {'probability': 'scattering_probability'}
     def apply(self, theory):
         if self.is2d:
 …
         Iq_calc = Iq_calc.reshape(self.nq, self.nq)
+        # CRUFT: don't need probability as a function anymore
         probability = self.probability() if callable(self.probability) else self.probability
         coverage = self.coverage

sasmodels/sasview_model.py

-                      ra8a1d48
+                      r9150036
 from . import core
 from . import custom
+from . import kernelcl
 from . import product
 from . import generate
 …
 from . import modelinfo
 from .details import make_kernel_args, dispersion_mesh
+from .kernelcl import reset_environment
 # pylint: disable=unused-import
 …
 #: has changed since we last reloaded.
 _CACHED_MODULE = {}  # type: Dict[str, "module"]
+def reset_environment():
+    # type: () -> None
+    """
+    Clear the compute engine context so that the GUI can change devices.
+    This removes all compiled kernels, even those that are active on fit
+    pages, but they will be restored the next time they are needed.
+    """
+    kernelcl.reset_environment()
+    for model in MODELS.values():
+        model._model = None
 def find_model(modelname):
 …
     def _calculate_Iq(self, qx, qy=None):
         if self._model is None:
+            self._model = core.build_model(self._model_info)
+            # Only need one copy of the compiled kernel regardless of how many
+            # times it is used, so store it in the class.  Also, to reset the
+            # compute engine, need to clear out all existing compiled kernels,
+            # which is much easier to do if we store them in the class.
+            self.__class__._model = core.build_model(self._model_info)
         if qy is not None:
             q_vectors = [np.asarray(qx), np.asarray(qy)]

README.rst

-                      re30d645
+                      r2a64722
 is available.
+Example
+Install
 -------
+The easiest way to use sasmodels is from `SasView <http://www.sasview.org/>`_.
+You can also install sasmodels as a standalone package in python. Use
+`miniconda <https://docs.conda.io/en/latest/miniconda.html>`_
+or `anaconda <https://www.anaconda.com/>`_
+to create a python environment with the sasmodels dependencies::
+    $ conda create -n sasmodels -c conda-forge numpy scipy matplotlib pyopencl
+The option ``-n sasmodels`` names the environment sasmodels, and the option
+``-c conda-forge`` selects the conda-forge package channel because pyopencl
+is not part of the base anaconda distribution.
+Activate the environment and install sasmodels::
+    $ conda activate sasmodels
+    (sasmodels) $ pip install sasmodels
+Install `bumps <https://github.com/bumps/bumps>`_ if you want to use it to fit
+your data::
+    (sasmodels) $ pip install bumps
+Usage
+-----
+Check that the works::
+    (sasmodels) $ python -m sasmodels.compare cylinder
+To show the orientation explorer::
+    (sasmodels) $ python -m sasmodels.jitter
+Documentation is available online as part of the SasView
+`fitting perspective <http://www.sasview.org/docs/index.html>`_
+as well as separate pages for
+`individual models <http://www.sasview.org/docs/user/sasgui/perspectives/fitting/models/index.html>`_.
+Programming details for sasmodels are available in the
+`developer documentation <http://www.sasview.org/docs/dev/dev.html>`_.
+Fitting Example
+---------------
 The example directory contains a radial+tangential data set for an oriented
 rod-like shape.
+The data is loaded by sas.dataloader from the sasview package, so sasview
+is needed to run the example.
+To load the example data, you will need the SAS data loader from the sasview
+package. This is not yet available on PyPI, so you will need a copy of the
+SasView source code to run it.  Create a directory somewhere to hold the
+sasview and sasmodels source code, which we will refer to as $SOURCE.
+To run the example, you need sasview, sasmodels and bumps.  Assuming these
+repositories are installed side by side, change to the sasmodels/example
+directory and enter::
+Use the following to install sasview, and the sasmodels examples::
+    PYTHONPATH=..:../../sasview/src ../../bumps/run.py fit.py \
+        cylinder --preview
+    (sasmodels) $ cd $SOURCE
+    (sasmodels) $ conda install git
+    (sasmodels) $ git clone https://github.com/sasview/sasview.git
+    (sasmodels) $ git clone https://github.com/sasview/sasmodels.git
+See bumps documentation for instructions on running the fit.  With the
+python packages installed, e.g., into a virtual environment, then the
+python path need not be set, and the command would be::
+Set the path to the sasview source on your python path within the sasmodels
+environment.  On Windows, this will be::
+    bumps fit.py cylinder --preview
+    (sasmodels)> set PYTHONPATH="$SOURCE\sasview\src"
+    (sasmodels)> cd $SOURCE/sasmodels/example
+    (sasmodels)> python -m bumps.cli fit.py cylinder --preview
+On Mac/Linux with the standard shell this will be::
+    (sasmodels) $ export PYTHONPATH="$SOURCE/sasview/src"
+    (sasmodels) $ cd $SOURCE/sasmodels/example
+    (sasmodels) $ bumps fit.py cylinder --preview
 The fit.py model accepts up to two arguments.  The first argument is the
 …
 both radial and tangential simultaneously, use the word "both".
+Notes
+-----
+cylinder.c + cylinder.py is the cylinder model with renamed variables and
+sld scaled by 1e6 so the numbers are nicer.  The model name is "cylinder"
+lamellar.py is an example of a single file model with embedded C code.
+See `bumps documentation <https://bumps.readthedocs.io/>`_ for detailed
+instructions on running the fit.
 |TravisStatus|_

doc/guide/gpu_setup.rst

-                      r63602b1
+                      r8b31efa
 Device Selection
 ================
+**OpenCL drivers**
 If you have multiple GPU devices you can tell the program which device to use.
 By default, the program looks for one GPU and one CPU device from available
 …
 was used to run the model.
+**If you don't want to use OpenCL, you can set** *SAS_OPENCL=None*
+**in your environment settings, and it will only use normal programs.**
+If you want to use one of the other devices, you can run the following
+If you want to use a specific driver and devices, you can run the following
 from the python console::
 …
 This will provide a menu of different OpenCL drivers available.
 When one is selected, it will say "set PYOPENCL_CTX=..."
+Use that value as the value of *SAS_OPENCL*.
+Use that value as the value of *SAS_OPENCL=driver:device*.
+To use the default OpenCL device (rather than CUDA or None),
+set *SAS_OPENCL=opencl*.
+In batch queues, you may need to set *XDG_CACHE_HOME=~/.cache*
+(Linux only) to a different directory, depending on how the filesystem
+is configured.  You should also set *SAS_DLL_PATH* for CPU-only modules.
+    -DSAS_MODELPATH=path sets directory containing custom models
+    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device
+    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only)
+    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler
+    -DSAS_OPENMP=1 turns on OpenMP for the DLLs
+    -DSAS_DLL_PATH=path sets the path to the compiled modules
+**CUDA drivers**
+If OpenCL drivers are not available on your system, but NVidia CUDA
+drivers are available, then set *SAS_OPENCL=cuda* or
+*SAS_OPENCL=cuda:n* for a particular device number *n*.  If no device
+number is specified, then the CUDA drivers looks for look for
+*CUDA_DEVICE=n* or a file ~/.cuda-device containing n for the device number.
+In batch queues, the SLURM command *sbatch --gres=gpu:1 ...* will set
+*CUDA_VISIBLE_DEVICES=n*, which ought to set the correct device
+number for *SAS_OPENCL=cuda*.  If not, then set
+*CUDA_DEVICE=$CUDA_VISIBLE_DEVICES* within the batch script.  You may
+need to set the CUDA cache directory to a folder accessible across the
+cluster with *PYCUDA_CACHE_DIR* (or *PYCUDA_DISABLE_CACHE* to disable
+caching), and you may need to set environment specific compiler flags
+with *PYCUDA_DEFAULT_NVCC_FLAGS*.  You should also set *SAS_DLL_PATH*
+for CPU-only modules.
+**No GPU support**
+If you don't want to use OpenCL or CUDA, you can set *SAS_OPENCL=None*
+in your environment settings, and it will only use normal programs.
+In batch queues, you may need to set *SAS_DLL_PATH* to a directory
+accessible on the compute node.
 Device Testing
 …
 *Document History*
 | 2017-09-27 Paul Kienzle
+| 2018-10-15 Paul Kienzle

doc/guide/scripting.rst

-                      rbd7630d
+                      r23df833
 python kernel.  Once the kernel is in hand, we can then marshal a set of
 parameters into a :class:`sasmodels.details.CallDetails` object and ship it to
+the kernel using the :func:`sansmodels.direct_model.call_kernel` function.  An
+example should help, *example/cylinder_eval.py*::
+    from numpy import logspace
+the kernel using the :func:`sansmodels.direct_model.call_kernel` function.  To
+accesses the underlying $<F(q)>$ and $<F^2(q)>$, use
+:func:`sasmodels.direct_model.call_Fq` instead.
+The following example should
+help, *example/cylinder_eval.py*::
+    from numpy import logspace, sqrt
     from matplotlib import pyplot as plt
     from sasmodels.core import load_model
     from sasmodels.direct_model import call_kernel
+    from sasmodels.direct_model import call_kernel, call_Fq
     model = load_model('cylinder')
     q = logspace(-3, -1, 200)
     kernel = model.make_kernel([q])
+    Iq = call_kernel(kernel, dict(radius=200.))
+    plt.loglog(q, Iq)
+    pars = {'radius': 200, 'radius_pd': 0.1, 'scale': 2}
+    Iq = call_kernel(kernel, pars)
+    F, Fsq, Reff, V, Vratio = call_Fq(kernel, pars)
+    plt.loglog(q, Iq, label='2 I(q)')
+    plt.loglog(q, F**2/V, label='<F(q)>^2/V')
+    plt.loglog(q, Fsq/V, label='<F^2(q)>/V')
+    plt.xlabel('q (1/A)')
+    plt.ylabel('I(q) (1/cm)')
+    plt.title('Cylinder with radius 200.')
+    plt.legend()
     plt.show()
+On windows, this can be called from the cmd prompt using sasview as::
+.. figure:: direct_call.png
+    Comparison between $I(q)$, $<F(q)>$ and $<F^2(q)>$ for cylinder model.
+This compares $I(q)$ with $<F(q)>$ and $<F^2(q)>$ for a cylinder
+with *radius=200 +/- 20* and *scale=2*. Note that *call_Fq* does not
+include scale and background, nor does it normalize by the average volume.
+The definition of $F = \rho V \hat F$ scaled by the contrast and
+volume, compared to the canonical cylinder $\hat F$, with $\hat F(0) = 1$.
+Integrating over polydispersity and orientation, the returned values are
+$\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F(q,r_o,\theta)\sin\theta$ and
+$\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F^2(q,r_o,\theta)\sin\theta$.
+On windows, this example can be called from the cmd prompt using sasview as
+as the python interpreter::
     SasViewCom example/cylinder_eval.py

example/cylinder_eval.py

-                      r2e66ef5
+                      r23df833
 """
 from numpy import logspace
+from numpy import logspace, sqrt
 from matplotlib import pyplot as plt
 from sasmodels.core import load_model
 from sasmodels.direct_model import call_kernel
+from sasmodels.direct_model import call_kernel, call_Fq
 model = load_model('cylinder')
 q = logspace(-3, -1, 200)
 kernel = model.make_kernel([q])
+Iq = call_kernel(kernel, dict(radius=200.))
+plt.loglog(q, Iq)
+pars = {'radius': 200, 'radius_pd': 0.1, 'scale': 2}
+Iq = call_kernel(kernel, pars)
+F, Fsq, Reff, V, Vratio = call_Fq(kernel, pars)
+plt.loglog(q, Iq, label='2 I(q)')
+plt.loglog(q, F**2/V, label='<F(q)>^2/V')
+plt.loglog(q, Fsq/V, label='<F^2(q)>/V')
 plt.xlabel('q (1/A)')
 plt.ylabel('I(q)')
+plt.ylabel('I(q) (1/cm)')
 plt.title('Cylinder with radius 200.')
+plt.legend()
 plt.show()

explore/precision.py

-                      rfba9ca0
+                      rcd28947
     return model_info
+# Hack to allow second parameter A in two parameter functions
+# Hack to allow second parameter A in the gammainc and gammaincc functions.
+# Create a 2-D variant of the precision test if we need to handle other two
+# parameter functions.
 A = 1
 def parse_extra_pars():
+    """
+    Parse the command line looking for the second parameter "A=..." for the
+    gammainc/gammaincc functions.
+    """
     global A
 …
+)
 add_function(
+    # Note: "a" is given as A=... on the command line via parse_extra_pars
     name="gammainc(x)",
     mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a),
 …
+)
 add_function(
+    # Note: "a" is given as A=... on the command line via parse_extra_pars
     name="gammaincc(x)",
     mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a),
 …
+)
 add_function(
     name="debye",
+    name="gauss_coil",
     mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
     np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4,
     ocl_function=make_ocl("""
     const double qsq = q*q;
+    if (qsq < 1.0) { // Pade approximation
+    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide)
+    // For single: use O(7) Taylor with 0.8 cutoff (7 mad)
+    if (qsq < 0.0) {
         const double x = qsq;
         if (0) { // 0.36 single
 …
             const double B1=3./8., B2=3./56., B3=1./336.;
             return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.);
         } else if (1) { // 1.0 for single, 0.25 for double
+        } else if (0) { // 1.0 for single, 0.25 for double
             // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
             const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
 …
                   /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
+        }
     } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double
+    } else if (qsq < 0.8) {
         const double x = qsq;
         const double C0 = +1.;

sasmodels/compare.py

-                      r610ef23
+                      rc1799d3
 from . import kerneldll
 from . import kernelcl
+from . import kernelcuda
 from .data import plot_theory, empty_data1D, empty_data2D, load_data
 from .direct_model import DirectModel, get_mesh
 …
     === environment variables ===
     -DSAS_MODELPATH=path sets directory containing custom models
     -DSAS_OPENCL=vendor:device|none sets the target OpenCL device
+    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device
     -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only)
     -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler
 …
     # If it is a list of choices, pick one at random with equal probability
-    # In practice, the model specific random generator will override.
     par = model_info.parameters[name]
     if len(par.limits) > 2:  # choice list
         return np.random.randint(len(par.limits))
+    if par.choices:  # choice list
+        return np.random.randint(len(par.choices))
     # If it is a fixed range, pick from it with equal probability.
 …
         set_integration_size(model_info, ngauss)
+    if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl():
+    if (dtype != "default" and not dtype.endswith('!')
+            and not (kernelcl.use_opencl() or kernelcuda.use_cuda())):
         raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR)
 …
         'rel_err'   : True,
         'explore'   : False,
         'use_demo'  : True,
+        'use_demo'  : False,
         'zero'      : False,
         'html'      : False,

sasmodels/core.py

-                      r2dcd6e7
+                      r07646b6
 from . import mixture
 from . import kernelpy
+from . import kernelcuda
 from . import kernelcl
 from . import kerneldll
 …
     numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform)
     source = generate.make_source(model_info)
     if platform == "dll":
         #print("building dll", numpy_dtype)
         return kerneldll.load_dll(source['dll'], model_info, numpy_dtype)
+    elif platform == "cuda":
+        return kernelcuda.GpuModel(source, model_info, numpy_dtype, fast=fast)
     else:
         #print("building ocl", numpy_dtype)
 …
     # type: (ModelInfo, str, str) -> (np.dtype, bool, str)
     """
     Interpret dtype string, returning np.dtype and fast flag.
+    Interpret dtype string, returning np.dtype, fast flag and platform.
     Possible types include 'half', 'single', 'double' and 'quad'.  If the
 …
     default for the model and platform.
+    Platform preference can be specfied ("ocl" vs "dll"), with the default
+    being OpenCL if it is availabe.  If the dtype name ends with '!' then
+    platform is forced to be DLL rather than OpenCL.
+    Platform preference can be specfied ("ocl", "cuda", "dll"), with the
+    default being OpenCL or CUDA if available, otherwise DLL.  If the dtype
+    name ends with '!' then platform is forced to be DLL rather than GPU.
+    The default platform is set by the environment variable SAS_OPENCL,
+    SAS_OPENCL=driver:device for OpenCL, SAS_OPENCL=cuda:device for CUDA
+    or SAS_OPENCL=none for DLL.
     This routine ignores the preferences within the model definition.  This
 …
     # Assign default platform, overriding ocl with dll if OpenCL is unavailable
     # If opencl=False OpenCL is switched off
     if platform is None:
         platform = "ocl"
-    if not kernelcl.use_opencl() or not model_info.opencl:
-        platform = "dll"
     # Check if type indicates dll regardless of which platform is given
 …
         platform = "dll"
         dtype = dtype[:-1]
+    # Make sure model allows opencl/gpu
+    if not model_info.opencl:
+        platform = "dll"
+    # Make sure opencl is available, or fallback to cuda then to dll
+    if platform == "ocl" and not kernelcl.use_opencl():
+        platform = "cuda" if kernelcuda.use_cuda() else "dll"
     # Convert special type names "half", "fast", and "quad"
 …
         dtype = "float16"
+    # Convert dtype string to numpy dtype.
+    # Convert dtype string to numpy dtype.  Use single precision for GPU
+    # if model allows it, otherwise use double precision.
     if dtype is None or dtype == "default":
         numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single
+        numpy_dtype = (generate.F32 if model_info.single and platform in ("ocl", "cuda")
                        else generate.F64)
     else:
         numpy_dtype = np.dtype(dtype)
     # Make sure that the type is supported by opencl, otherwise use dll
+    # Make sure that the type is supported by GPU, otherwise use dll
     if platform == "ocl":
         env = kernelcl.environment()
+        if not env.has_type(numpy_dtype):
+            platform = "dll"
+            if dtype is None:
+                numpy_dtype = generate.F64
+    elif platform == "cuda":
+        env = kernelcuda.environment()
+    else:
+        env = None
+    if env is not None and not env.has_type(numpy_dtype):
+        platform = "dll"
+        if dtype is None:
+            numpy_dtype = generate.F64
     return numpy_dtype, fast, platform
 …
     model = load_model("cylinder@hardsphere*sphere")
     actual = [p.name for p in model.info.parameters.kernel_parameters]
+    target = ("sld sld_solvent radius length theta phi volfraction"
+    target = ("sld sld_solvent radius length theta phi"
+              " radius_effective volfraction "
+              " structure_factor_mode radius_effective_mode"
               " A_sld A_sld_solvent A_radius").split()
     assert target == actual, "%s != %s"%(target, actual)

sasmodels/data.py

rbd7630d	rc88f983
337	337	else:
338	338	dq = resolution * q
339
340	339	data = Data1D(q, Iq, dx=dq, dy=dIq)
341	340	data.filename = "fake data"

sasmodels/generate.py

-                      r6e45516
+                      ra8a1d48
     dimension, or 1.0 if no volume normalization is required.
+    *ER(p1, p2, ...)* returns the effective radius of the form with
+    particular dimensions.
+    *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms.
+    *shell_volume(p1, p2, ...)* returns the volume of the shell for forms
+    which are hollow.
+    *effective_radius(mode, p1, p2, ...)* returns the effective radius of
+    the form with particular dimensions.  Mode determines the type of
+    effective radius returned, with mode=1 for equivalent volume.
     #define INVALID(v) (expr)  returns False if v.parameter is invalid
 …
 background.  This will work correctly even when polydispersity is off.
-*ER* and *VR* are python functions which operate on parameter vectors.
-The constructor code will generate the necessary vectors for computing
-them with the desired polydispersity.
 The kernel module must set variables defining the kernel meta data:
 …
     create the OpenCL kernel functions.  The files defining the functions
     need to be listed before the files which use the functions.
-    *ER* is a python function defining the effective radius.  If it is
-    not present, the effective radius is 0.
-    *VR* is a python function defining the volume ratio.  If it is not
-    present, the volume ratio is 1.
     *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing
 …
+# type in IQXY pattern could be single, float, double, long double, ...
+_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]",
+_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]",
                            flags=re.MULTILINE)
 def find_xy_mode(source):
 …
     return 'qa'
+# Note: not presently used.  Need to know whether Fq is available before
+# trying to compile the source.  Leave the code here in case we decide that
+# define have_Fq for each form factor is too tedious and error prone.
+_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE)
+def contains_Fq(source):
+    # type: (List[str]) -> bool
+    """
+    Return True if C source defines "void Fq(".
+    """
+    for code in source:
+        if _FQ_PATTERN.search(code) is not None:
+            return True
+    return False
+_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE)
+def contains_shell_volume(source):
+    # type: (List[str]) -> bool
+    """
+    Return True if C source defines "double shell_volume(".
+    """
+    for code in source:
+        if _SHELL_VOLUME_PATTERN.search(code) is not None:
+            return True
+    return False
 def _add_source(source, code, path, lineno=1):
 …
     # dispersion.  Need to be careful that necessary parameters are available
     # for computing volume even if we allow non-disperse volume parameters.
     partable = model_info.parameters
 …
     for path, code in user_code:
         _add_source(source, code, path)
     if model_info.c_code:
         _add_source(source, model_info.c_code, model_info.filename,
 …
     q, qx, qy, qab, qa, qb, qc \
         = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
     # Generate form_volume function, etc. from body only
     if isinstance(model_info.form_volume, str):
         pars = partable.form_volume_parameters
         source.append(_gen_fn(model_info, 'form_volume', pars))
+    if isinstance(model_info.shell_volume, str):
+        pars = partable.form_volume_parameters
+        source.append(_gen_fn(model_info, 'shell_volume', pars))
     if isinstance(model_info.Iq, str):
         pars = [q] + partable.iq_parameters
 …
         pars = [qa, qb, qc] + partable.iq_parameters
         source.append(_gen_fn(model_info, 'Iqabc', pars))
+    # Check for shell_volume in source
+    is_hollow = contains_shell_volume(source)
     # What kind of 2D model do we need?  Is it consistent with the parameters?
 …
     source.append("\\\n".join(p.as_definition()
                               for p in partable.kernel_parameters))
     # Define the function calls
+    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0"
     if partable.form_volume_parameters:
         refs = _call_pars("_v.", partable.form_volume_parameters)
+        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs))
+        if is_hollow:
+            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2)
+        else:
+            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs))
+        if model_info.effective_radius_type:
+            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %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"
+        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)"
     source.append(call_volume)
+    source.append(call_effective_radius)
     model_refs = _call_pars("_v.", partable.iq_parameters)
+    pars = ",".join(["_q"] + model_refs)
+    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
+    if model_info.have_Fq:
+        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
+        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars
+        clear_iq = "#undef CALL_FQ"
+    else:
+        pars = ",".join(["_q"] + model_refs)
+        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
+        clear_iq = "#undef CALL_IQ"
     if xy_mode == 'qabc':
         pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
 …
         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_AC"
     elif xy_mode == 'qa':
+    elif xy_mode == 'qa' and not model_info.have_Fq:
         pars = ",".join(["_qa"] + model_refs)
         call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_A"
+    elif xy_mode == 'qa' and model_info.have_Fq:
+        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs)
+        # Note: uses rare C construction (expr1, expr2) which computes
+        # expr1 then expr2 and evaluates to expr2.  This allows us to
+        # leave it looking like a function even though it is returning
+        # its values by reference.
+        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars
+        clear_iqxy = "#undef CALL_FQ_A"
     elif xy_mode == 'qxy':
         orientation_refs = _call_pars("_v.", partable.orientation_parameters)
 …
     magpars = [k-2 for k, p in enumerate(partable.call_parameters)
                if p.type == 'sld']
     # Fill in definitions for numbers of parameters
     source.append("#define MAX_PD %s"%partable.max_pd)
 …
     source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
     source.append("#define PROJECTION %d"%PROJECTION)
     # TODO: allow mixed python/opencl kernels?
     ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
+    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
+    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
+    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
     result = {
         'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
         'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
+    }
     return result
 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name):
+def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name):
     # type: ([str,str], str, str, str) -> List[str]
     code = kernel[0]
 …
         '#line 1 "%s Iq"' % path,
         code,
         "#undef CALL_IQ",
+        clear_iq,
         "#undef KERNEL_NAME",
+        ]
 …
         pars = model_info.parameters.kernel_parameters
     else:
+        pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters
+        pars = (model_info.parameters.common_parameters
+                + model_info.parameters.kernel_parameters)
     partable = make_partable(pars)
     subst = dict(id=model_info.id.replace('_', '-'),

sasmodels/jitter.py

-                      r1198f90
+                      r7d97437
     pass
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
 …
         pass
+    axcolor = 'lightgoldenrodyellow'
+    # CRUFT: use axisbg instead of facecolor for matplotlib<2
+    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg'
+    props = {facecolor_prop: 'lightgoldenrodyellow'}
     ## add control widgets to plot
     axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)
     axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor)
     axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor)
+    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props)
+    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props)
+    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props)
     stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta)
     sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi)
     spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi)
     axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)
     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor)
     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor)
+    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props)
+    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props)
+    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props)
     # Note: using ridiculous definition of rectangle distribution, whose width
     # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep

sasmodels/kernel.py

-                      r2d81cfe
+                      rcd28947
     dtype = None  # type: np.dtype
     def __call__(self, call_details, values, cutoff, magnetic):
+    def Iq(self, call_details, values, cutoff, magnetic):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        raise NotImplementedError("need to implement __call__")
+        r"""
+        Returns I(q) from the polydisperse average scattering.
+        .. math::
+            I(q) = \text{scale} \cdot P(q) + \text{background}
+        With the correct choice of model and contrast, setting *scale* to
+        the volume fraction $V_f$ of particles should match the measured
+        absolute scattering.  Some models (e.g., vesicle) have volume fraction
+        built into the model, and do not need an additional scale.
+        """
+        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic,
+                              effective_radius_type=0)
+        combined_scale = values[0]/shell_volume
+        background = values[1]
+        return combined_scale*F2 + background
+    __call__ = Iq
+    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        r"""
+        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and
+        form:shell volume ratio.  The <F(q)> term may be None if the
+        form factor does not support direct computation of $F(q)$
+        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations,
+        .. math::
+            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background}
+        For the beta approximation, this becomes
+        .. math::
+            I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background}
+                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background}
+        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape
+        and orientation, with each configuration $x_k$ having form factor
+        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is:
+        .. math::
+            P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k}
+        The form factor itself is scaled by volume and contrast to compute the
+        total scattering.  This is then squared, and the volume weighted
+        F^2 is then normalized by volume F.  For a given density, the number
+        of scattering centers is assumed to scale linearly with volume.  Later
+        scaling the resulting $P(q)$ by the volume fraction of particles
+        gives the total scattering on an absolute scale. Most models
+        incorporate the volume fraction into the overall scale parameter.  An
+        exception is vesicle, which includes the volume fraction parameter in
+        the model itself, scaling $F$ by $\surd V_f$ so that the math for
+        the beta approximation works out.
+        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make
+        sure that the polydisperisity distributions normalize to one.  In
+        particular, any distibution values $x_k$ outside the valid domain
+        of $F$ will not be included, and the distribution will be implicitly
+        truncated.  This is controlled by the parameter limits defined in the
+        model (which truncate the distribution before calling the kernel) as
+        well as any region excluded using the *INVALID* macro defined within
+        the model itself.
+        The volume used in the polydispersity calculation is the form volume
+        for solid objects or the shell volume for hollow objects.  Shell
+        volume should be used within $F$ so that the normalizing scale
+        represents the volume fraction of the shell rather than the entire
+        form.  This corresponds to the volume fraction of shell-forming
+        material added to the solvent.
+        The calculation of $S$ requires the effective radius and the
+        volume fraction of the particles.  The model can have several
+        different ways to compute effective radius, with the
+        *effective_radius_type* parameter used to select amongst them.  The
+        volume fraction of particles should be determined from the total
+        volume fraction of the form, not just the shell volume fraction.
+        This makes a difference for hollow shapes, which need to scale
+        the volume fraction by the returned volume ratio when computing $S$.
+        For solid objects, the shell volume is set to the form volume so
+        this scale factor evaluates to one and so can be used for both
+        hollow and solid shapes.
+        """
+        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type)
+        #print("returned",self.q_input.q, self.result)
+        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
+        total_weight = self.result[nout*self.q_input.nq + 0]
+        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it
+        # is okay to test directly against zero.  If weight is zero then I(q),
+        # etc. must also be zero.
+        if total_weight == 0.:
+            total_weight = 1.
+        # Note: shell_volume == form_volume for solid objects
+        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight
+        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight
+        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight
+        if shell_volume == 0.:
+            shell_volume = 1.
+        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None
+        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight
+        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume
     def release(self):
         # type: () -> None
         pass
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        """
+        Call the kernel.  Subclasses defining kernels for particular execution
+        engines need to provide an implementation for this.
+        """
+        raise NotImplementedError()

sasmodels/kernel_header.c

-                      r108e70e
+                      r07646b6
 #ifdef __OPENCL_VERSION__
 # define USE_OPENCL
+#elif defined(__CUDACC__)
+# define USE_CUDA
 #elif defined(_OPENMP)
 # define USE_OPENMP
 #endif
+// Use SAS_DOUBLE to force the use of double even for float kernels
+#define SAS_DOUBLE dou ## ble
 // If opencl is not available, then we are compiling a C function
 // Note: if using a C++ compiler, then define kernel as extern "C"
 #ifdef USE_OPENCL
+   #define USE_GPU
+   #define pglobal global
+   #define pconstant constant
    typedef int int32_t;
+#  if defined(USE_SINCOS)
+#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+#  else
+#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif
+   #if defined(USE_SINCOS)
+   #  define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+   #else
+   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+   #endif
    // Intel CPU on Mac gives strange values for erf(); on the verified
    // platforms (intel, nvidia, amd), the cephes erf() is significantly
 …
    #  define erfcf erfc
    #endif
+#else // !USE_OPENCL
+// Use SAS_DOUBLE to force the use of double even for float kernels
+#  define SAS_DOUBLE dou ## ble
+#  ifdef __cplusplus
+#elif defined(USE_CUDA)
+   #define USE_GPU
+   #define local __shared__
+   #define pglobal
+   #define constant __constant__
+   #define pconstant const
+   #define kernel extern "C" __global__
+   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+   // OpenCL pown(a,b) = C99 pow(a,b), b integer
+   #define powr(a,b) pow(a,b)
+   #define pown(a,b) pow(a,b)
+   //typedef int int32_t;
+   #if defined(USE_SINCOS)
+   #  define SINCOS(angle,svar,cvar) sincos(angle,&svar,&cvar)
+   #else
+   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+   #endif
+#else // !USE_OPENCL && !USE_CUDA
+   #define local
+   #define pglobal
+   #define constant const
+   #define pconstant const
+   #ifdef __cplusplus
       #include <cstdio>
       #include <cmath>
 …
      #endif
      inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
 #  else // !__cplusplus
+   #else // !__cplusplus
      #include <inttypes.h>  // C99 guarantees that int32_t types is here
      #include <stdio.h>
 …
          #define NEED_EXPM1
          #define NEED_TGAMMA
+         #define NEED_CBRT
          // expf missing from windows?
          #define expf exp
 …
      #define kernel
      #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif  // !__cplusplus
+#  define global
+#  define local
+#  define constant const
+// OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+// OpenCL pown(a,b) = C99 pow(a,b), b integer
+#  define powr(a,b) pow(a,b)
+#  define pown(a,b) pow(a,b)
+   #endif  // !__cplusplus
+   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0
+   // OpenCL pown(a,b) = C99 pow(a,b), b integer
+   #define powr(a,b) pow(a,b)
+   #define pown(a,b) pow(a,b)
 #endif // !USE_OPENCL
+#if defined(NEED_CBRT)
+   #define cbrt(_x) pow(_x, 0.33333333333333333333333)
+#endif
 #if defined(NEED_EXPM1)

sasmodels/kernel_iq.c

-                      r70530778
+                      r12f4c19
 //  MAGNETIC_PARS : a comma-separated list of indices to the sld
 //      parameters in the parameter table.
+//  CALL_VOLUME(table) : call the form volume function
+//  CALL_VOLUME(form, shell, table) : assign form and shell values
+//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function
 //  CALL_IQ(q, table) : call the Iq function for 1D calcs.
 //  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data.
+//  CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs.
+//  CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data.
 //  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes
 //  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes
 …
 //  PROJECTION : equirectangular=1, sinusoidal=2
 //      see explore/jitter.py for definitions.
 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
 …
 #endif // _QABC_SECTION
 // ==================== KERNEL CODE ========================
 kernel
 void KERNEL_NAME(
+    int32_t nq,                 // number of q values
+    const int32_t pd_start,     // where we are in the dispersity loop
+    const int32_t pd_stop,      // where we are stopping in the dispersity loop
+    global const ProblemDetails *details,
+    global const double *values,
+    global const double *q, // nq q values, with padding to boundary
+    global double *result,  // nq+1 return values, again with padding
+    const double cutoff     // cutoff in the dispersity weight product
+    int32_t nq,                   // number of q values
+    const int32_t pd_start,       // where we are in the dispersity loop
+    const int32_t pd_stop,        // where we are stopping in the dispersity loop
+    pglobal const ProblemDetails *details,
+    pglobal const double *values, // parameter values and distributions
+    pglobal const double *q,      // nq q values, with padding to boundary
+    pglobal double *result,       // nq+1 return values, again with padding
+    const double cutoff,          // cutoff in the dispersity weight product
+    int32_t effective_radius_type // which effective radius to compute
+    )
+{
 #ifdef USE_OPENCL
+#if defined(USE_GPU)
   // who we are and what element we are working with
+  #if defined(USE_OPENCL)
   const int q_index = get_global_id(0);
+  #else // USE_CUDA
+  const int q_index = threadIdx.x + blockIdx.x * blockDim.x;
+  #endif
   if (q_index >= nq) return;
 #else
 …
   //
   // The code differs slightly between opencl and dll since opencl is only
   // seeing one q value (stored in the variable "this_result") while the dll
+  // seeing one q value (stored in the variable "this_F2") while the dll
   // version must loop over all q.
+  #ifdef USE_OPENCL
+    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
+  #else // !USE_OPENCL
+    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);
+  #if defined(CALL_FQ)
+    double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
+    double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+    double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]);
+    double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]);
+  #else
+    double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
+    double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]);
+    double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]);
+    double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]);
+  #endif
+  #if defined(USE_GPU)
+    #if defined(CALL_FQ)
+      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]);
+      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]);
+    #else
+      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]);
+    #endif
+  #else // !USE_GPU
     if (pd_start == 0) {
       #ifdef USE_OPENMP
       #pragma omp parallel for
       #endif
+      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
+      #if defined(CALL_FQ)
+          // 2*nq for F^2,F pairs
+          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0;
+      #else
+          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
+      #endif
+    }
     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
 #endif // !USE_OPENCL
+#endif // !USE_GPU
 …
   const int n4 = pd_length[4];
   const int p4 = pd_par[4];
   global const double *v4 = pd_value + pd_offset[4];
   global const double *w4 = pd_weight + pd_offset[4];
+  pglobal const double *v4 = pd_value + pd_offset[4];
+  pglobal const double *w4 = pd_weight + pd_offset[4];
   int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start
 …
           FETCH_Q         // set qx,qy from the q input vector
           APPLY_ROTATION  // convert qx,qy to qa,qb,qc
           CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
+          CALL_KERNEL     // F2 = Iqxy(qa, qb, qc, p1, p2, ...)
       ++step;  // increment counter representing position in dispersity mesh
 …
 // inner loop and defines the macros that use them.
+#if defined(CALL_IQ)
+  // unoriented 1D
+#if defined(CALL_FQ) // COMPUTE_F1_F2 is true
+  // unoriented 1D returning <F> and <F^2>
+  // Note that F1 and F2 are returned from CALL_FQ by reference, and the
+  // user of the CALL_KERNEL macro below is assuming that F1 and F2 are defined.
+  double qk;
+  double F1, F2;
+  #define FETCH_Q() do { qk = q[q_index]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table)
+#elif defined(CALL_FQ_A)
+  // unoriented 2D return <F> and <F^2>
+  // Note that the CALL_FQ_A macro is computing _F1_slot and _F2_slot by
+  // reference then returning _F2_slot.  We are calling them _F1_slot and
+  // _F2_slot here so they don't conflict with _F1 and _F2 in the macro
+  // expansion, or with the use of F2 = CALL_KERNEL() when it is used below.
+  double qx, qy;
+  double _F1_slot, _F2_slot;
+  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0)
+  #define BUILD_ROTATION() do {} while(0)
+  #define APPLY_ROTATION() do {} while(0)
+  #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),_F1_slot,_F2_slot,local_values.table)
+#elif defined(CALL_IQ)
+  // unoriented 1D return <F^2>
   double qk;
   #define FETCH_Q() do { qk = q[q_index]; } while (0)
   #define BUILD_ROTATION() do {} while(0)
   #define APPLY_ROTATION() do {} while(0)
   #define CALL_KERNEL() CALL_IQ(qk, local_values.table)
+  #define CALL_KERNEL() CALL_IQ(qk,local_values.table)
 #elif defined(CALL_IQ_A)
 …
   #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc)
   #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table)
 #elif defined(CALL_IQ_XY)
   // direct call to qx,qy calculator
 …
 // logic in the IQ_AC and IQ_ABC branches.  This will become more important
 // if we implement more projections, or more complicated projections.
+#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation
+#if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A)
+  // no orientation
   #define APPLY_PROJECTION() const double weight=weight0
 #elif defined(CALL_IQ_XY) // pass orientation to the model
 …
   const int n##_LOOP = details->pd_length[_LOOP]; \
   const int p##_LOOP = details->pd_par[_LOOP]; \
   global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
   global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
+  pglobal const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
+  pglobal const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
   int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
 …
 // Pointers to the start of the dispersity and weight vectors, if needed.
 #if MAX_PD>0
   global const double *pd_value = values + NUM_VALUES;
   global const double *pd_weight = pd_value + details->num_weights;
+  pglobal const double *pd_value = values + NUM_VALUES;
+  pglobal const double *pd_weight = pd_value + details->num_weights;
 #endif
 …
     // Note: weight==0 must always be excluded
     if (weight > cutoff) {
+      pd_norm += weight * CALL_VOLUME(local_values.table);
+      double form, shell;
+      CALL_VOLUME(form, shell, local_values.table);
+      weight_norm += weight;
+      weighted_form += weight * form;
+      weighted_shell += weight * shell;
+      if (effective_radius_type != 0) {
+        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table);
+      }
       BUILD_ROTATION();
 #ifndef USE_OPENCL
+#if !defined(USE_GPU)
       // DLL needs to explicitly loop over the q values.
       #ifdef USE_OPENMP
 …
       #endif
       for (q_index=0; q_index<nq; q_index++)
 #endif // !USE_OPENCL
+#endif // !USE_GPU
+      {
 …
         #if defined(MAGNETIC) && NUM_MAGNETIC > 0
           // Compute the scattering from the magnetic cross sections.
           double scattering = 0.0;
+          double F2 = 0.0;
           const double qsq = qx*qx + qy*qy;
           if (qsq > 1.e-16) {
 …
 //  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz);
+                }
                 scattering += xs_weight * CALL_KERNEL();
+                F2 += xs_weight * CALL_KERNEL();
+              }
+            }
+          }
         #else  // !MAGNETIC
+          const double scattering = CALL_KERNEL();
+          #if defined(CALL_FQ)
+            CALL_KERNEL(); // sets F1 and F2 by reference
+          #else
+            const double F2 = CALL_KERNEL();
+          #endif
         #endif // !MAGNETIC
+//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0);
+        #ifdef USE_OPENCL
+          this_result += weight * scattering;
+//printf("q_index:%d %g %g %g %g\n", q_index, F2, weight0);
+        #if defined(USE_GPU)
+          #if defined(CALL_FQ)
+            this_F2 += weight * F2;
+            this_F1 += weight * F1;
+          #else
+            this_F2 += weight * F2;
+          #endif
         #else // !USE_OPENCL
+          result[q_index] += weight * scattering;
+          #if defined(CALL_FQ)
+            result[2*q_index+0] += weight * F2;
+            result[2*q_index+1] += weight * F1;
+          #else
+            result[q_index] += weight * F2;
+          #endif
         #endif // !USE_OPENCL
+      }
+    }
+  }
 // close nested loops
 ++step;
 …
 #endif
+// Remember the current result and the updated norm.
+#ifdef USE_OPENCL
+  result[q_index] = this_result;
+  if (q_index == 0) result[nq] = pd_norm;
+//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm);
+#else // !USE_OPENCL
+  result[nq] = pd_norm;
+//printf("res: %g/%g\n", result[0], pd_norm);
+#endif // !USE_OPENCL
+// Remember the results and the updated norm.
+#if defined(USE_GPU)
+  #if defined(CALL_FQ)
+  result[2*q_index+0] = this_F2;
+  result[2*q_index+1] = this_F1;
+  #else
+  result[q_index] = this_F2;
+  #endif
+  if (q_index == 0)
+#endif
+  {
+#if defined(CALL_FQ)
+    result[2*nq] = weight_norm;
+    result[2*nq+1] = weighted_form;
+    result[2*nq+2] = weighted_shell;
+    result[2*nq+3] = weighted_radius;
+#else
+    result[nq] = weight_norm;
+    result[nq+1] = weighted_form;
+    result[nq+2] = weighted_shell;
+    result[nq+3] = weighted_radius;
+#endif
+  }
 // ** clear the macros in preparation for the next kernel **

sasmodels/kernelcl.py

-                      rd86f0fc
+                      r8064d5e
 """
 GPU driver for C kernels
+TODO: docs are out of date
 There should be a single GPU environment running on the system.  This
 …
 import time
+try:
+    from time import perf_counter as clock
+except ImportError: # CRUFT: python < 3.3
+    import sys
+    if sys.platform.count("darwin") > 0:
+        from time import time as clock
+    else:
+        from time import clock
 import numpy as np  # type: ignore
+# Attempt to setup opencl. This may fail if the opencl package is not
+# Attempt to setup opencl. This may fail if the pyopencl package is not
 # installed or if it is installed but there are no devices available.
 try:
 …
 from . import generate
+from .generate import F32, F64
 from .kernel import KernelModel, Kernel
 …
 def use_opencl():
+    return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none"
+    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower()
+    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda")
 ENV = None
 …
     Return true if device supports the requested precision.
     """
     if dtype == generate.F32:
+    if dtype == F32:
         return True
     elif dtype == generate.F64:
+    elif dtype == F64:
         return "cl_khr_fp64" in device.extensions
-    elif dtype == generate.F16:
-        return "cl_khr_fp16" in device.extensions
     else:
+        # Not supporting F16 type since it isn't accurate enough
         return False
 …
         cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
         queue.device)
-def _stretch_input(vector, dtype, extra=1e-3, boundary=32):
-    # type: (np.ndarray, np.dtype, float, int) -> np.ndarray
-    """
-    Stretch an input vector to the correct boundary.
-    Performance on the kernels can drop by a factor of two or more if the
-    number of values to compute does not fall on a nice power of two
-    boundary.   The trailing additional vector elements are given a
-    value of *extra*, and so f(*extra*) will be computed for each of
-    them.  The returned array will thus be a subset of the computed array.
-    *boundary* should be a power of 2 which is at least 32 for good
-    performance on current platforms (as of Jan 2015).  It should
-    probably be the max of get_warp(kernel,queue) and
-    device.min_data_type_align_size//4.
-    """
-    remainder = vector.size % boundary
-    if remainder != 0:
-        size = vector.size + (boundary - remainder)
-        vector = np.hstack((vector, [extra] * (size - vector.size)))
-    return np.ascontiguousarray(vector, dtype=dtype)
 def compile_model(context, source, dtype, fast=False):
 …
     """
     GPU context, with possibly many devices, and one queue per device.
+    Because the environment can be reset during a live program (e.g., if the
+    user changes the active GPU device in the GUI), everything associated
+    with the device context must be cached in the environment and recreated
+    if the environment changes.  The *cache* attribute is a simple dictionary
+    which holds keys and references to objects, such as compiled kernels and
+    allocated buffers.  The running program should check in the cache for
+    long lived objects and create them if they are not there.  The program
+    should not hold onto cached objects, but instead only keep them active
+    for the duration of a function call.  When the environment is destroyed
+    then the *release* method for each active cache item is called before
+    the environment is freed.  This means that each cl buffer should be
+    in its own cache entry.
     """
     def __init__(self):
         # type: () -> None
         # find gpu context
+        #self.context = cl.create_some_context()
+        self.context = None
+        if 'SAS_OPENCL' in os.environ:
+            #Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context
+            os.environ["PYOPENCL_CTX"] = os.environ["SAS_OPENCL"]
+        if 'PYOPENCL_CTX' in os.environ:
+            self._create_some_context()
+        if not self.context:
+            self.context = _get_default_context()
+        context_list = _create_some_context()
+        # Find a context for F32 and for F64 (maybe the same one).
+        # F16 isn't good enough.
+        self.context = {}
+        for dtype in (F32, F64):
+            for context in context_list:
+                if has_type(context.devices[0], dtype):
+                    self.context[dtype] = context
+                    break
+            else:
+                self.context[dtype] = None
+        # Build a queue for each context
+        self.queue = {}
+        context = self.context[F32]
+        self.queue[F32] = cl.CommandQueue(context, context.devices[0])
+        if self.context[F64] == self.context[F32]:
+            self.queue[F64] = self.queue[F32]
+        else:
+            context = self.context[F64]
+            self.queue[F64] = cl.CommandQueue(context, context.devices[0])
         # Byte boundary for data alignment
         #self.data_boundary = max(d.min_data_type_align_size
         #                         for d in self.context.devices)
+        self.queues = [cl.CommandQueue(context, context.devices[0])
                        for context in self.context]
+        #self.data_boundary = max(context.devices[0].min_data_type_align_size
+        #                         for context in self.context.values())
+        # Cache for compiled programs, and for items in context
         self.compiled = {}
+        self.cache = {}
     def has_type(self, dtype):
 …
         Return True if all devices support a given type.
         """
+        return any(has_type(d, dtype)
+                   for context in self.context
+                   for d in context.devices)
+    def get_queue(self, dtype):
+        # type: (np.dtype) -> cl.CommandQueue
+        """
+        Return a command queue for the kernels of type dtype.
+        """
+        for context, queue in zip(self.context, self.queues):
+            if all(has_type(d, dtype) for d in context.devices):
+                return queue
+    def get_context(self, dtype):
+        # type: (np.dtype) -> cl.Context
+        """
+        Return a OpenCL context for the kernels of type dtype.
+        """
+        for context in self.context:
+            if all(has_type(d, dtype) for d in context.devices):
+                return context
+    def _create_some_context(self):
+        # type: () -> cl.Context
+        """
+        Protected call to cl.create_some_context without interactivity.  Use
+        this if SAS_OPENCL is set in the environment.  Sets the *context*
+        attribute.
+        """
+        try:
+            self.context = [cl.create_some_context(interactive=False)]
+        except Exception as exc:
+            warnings.warn(str(exc))
+            warnings.warn("pyopencl.create_some_context() failed")
+            warnings.warn("the environment variable 'SAS_OPENCL' might not be set correctly")
+        return self.context.get(dtype, None) is not None
     def compile_program(self, name, source, dtype, fast, timestamp):
 …
             del self.compiled[key]
         if key not in self.compiled:
             context = self.get_context(dtype)
+            context = self.context[dtype]
             logging.info("building %s for OpenCL %s", key,
                          context.devices[0].name.strip())
             program = compile_model(self.get_context(dtype),
+            program = compile_model(self.context[dtype],
                                     str(source), dtype, fast)
             self.compiled[key] = (program, timestamp)
         return program
+    def free_buffer(self, key):
+        if key in self.cache:
+            self.cache[key].release()
+            del self.cache[key]
+    def __del__(self):
+        for v in self.cache.values():
+            release = getattr(v, 'release', lambda: None)
+            release()
+        self.cache = {}
+_CURRENT_ID = 0
+def unique_id():
+    global _CURRENT_ID
+    _CURRENT_ID += 1
+    return _CURRENT_ID
+def _create_some_context():
+    # type: () -> cl.Context
+    """
+    Protected call to cl.create_some_context without interactivity.
+    Uses SAS_OPENCL or PYOPENCL_CTX if they are set in the environment,
+    otherwise scans for the most appropriate device using
+    :func:`_get_default_context`.  Ignore *SAS_OPENCL=OpenCL*, which
+    indicates that an OpenCL device should be used without specifying
+    which one (and not a CUDA device, or no GPU).
+    """
+    # Assume we do not get here if SAS_OPENCL is None or CUDA
+    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl')
+    if sas_opencl.lower() != 'opencl':
+        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context
+        os.environ["PYOPENCL_CTX"] = sas_opencl
+    if 'PYOPENCL_CTX' in os.environ:
+        try:
+            return [cl.create_some_context(interactive=False)]
+        except Exception as exc:
+            warnings.warn(str(exc))
+            warnings.warn("pyopencl.create_some_context() failed")
+            warnings.warn("the environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set correctly")
+    return _get_default_context()
 def _get_default_context():
 …
         self.dtype = dtype
         self.fast = fast
         self.program = None # delay program creation
         self._kernels = None
+        self.timestamp = generate.ocl_timestamp(self.info)
+        self._cache_key = unique_id()
     def __getstate__(self):
 …
         # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None
         self.info, self.source, self.dtype, self.fast = state
-        self.program = None
     def make_kernel(self, q_vectors):
         # type: (List[np.ndarray]) -> "GpuKernel"
+        if self.program is None:
+            compile_program = environment().compile_program
+            timestamp = generate.ocl_timestamp(self.info)
+            self.program = compile_program(
+        return GpuKernel(self, q_vectors)
+    @property
+    def Iq(self):
+        return self._fetch_kernel('Iq')
+    def fetch_kernel(self, name):
+        # type: (str) -> cl.Kernel
+        """
+        Fetch the kernel from the environment by name, compiling it if it
+        does not already exist.
+        """
+        gpu = environment()
+        key = self._cache_key
+        if key not in gpu.cache:
+            program = gpu.compile_program(
                 self.info.name,
                 self.source['opencl'],
                 self.dtype,
                 self.fast,
                 timestamp)
+                self.timestamp)
             variants = ['Iq', 'Iqxy', 'Imagnetic']
             names = [generate.kernel_name(self.info, k) for k in variants]
             kernels = [getattr(self.program, k) for k in names]
             self._kernels = dict((k, v) for k, v in zip(variants, kernels))
         is_2d = len(q_vectors) == 2
         if is_2d:
             kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']]
+            kernels = [getattr(program, k) for k in names]
+            data = dict((k, v) for k, v in zip(variants, kernels))
+            # keep a handle to program so GC doesn't collect
+            data['program'] = program
+            gpu.cache[key] = data
         else:
+            kernel = [self._kernels['Iq']]*2
+        return GpuKernel(kernel, self.dtype, self.info, q_vectors)
+    def release(self):
+        # type: () -> None
+        """
+        Free the resources associated with the model.
+        """
+        if self.program is not None:
+            self.program = None
+    def __del__(self):
+        # type: () -> None
+        self.release()
+            data = gpu.cache[key]
+        return data[name]
 # TODO: check that we don't need a destructor for buffers which go out of scope
 …
         # type: (List[np.ndarray], np.dtype) -> None
         # TODO: do we ever need double precision q?
-        env = environment()
         self.nq = q_vectors[0].size
         self.dtype = np.dtype(dtype)
 …
         # architectures tested so far.
         if self.is_2d:
+            # Note: 16 rather than 15 because result is 1 longer than input.
+            width = ((self.nq+16)//16)*16
+            width = ((self.nq+15)//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: 32 rather than 31 because result is 1 longer than input.
+            width = ((self.nq+32)//32)*32
+            width = ((self.nq+31)//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)
+        #print("creating inputs of size", self.global_size)
+        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                             hostbuf=self.q)
+        self._cache_key = unique_id()
+    @property
+    def q_b(self):
+        """Lazy creation of q buffer so it can survive context reset"""
+        env = environment()
+        key = self._cache_key
+        if key not in env.cache:
+            context = env.context[self.dtype]
+            #print("creating inputs of size", self.global_size)
+            buffer = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                               hostbuf=self.q)
+            env.cache[key] = buffer
+        return env.cache[key]
     def release(self):
         # type: () -> None
         """
+        Free the memory.
+        """
+        if self.q_b is not None:
+            self.q_b.release()
+            self.q_b = None
+        Free the buffer associated with the q value
+        """
+        environment().free_buffer(id(self))
     def __del__(self):
 …
     Callable SAS kernel.
     *kernel* is the GpuKernel object to call
     *model_info* is the module information
     *q_vectors* is the q vectors at which the kernel should be evaluated
+    *model* is the GpuModel object to call
+    The following attributes are defined:
+    *info* is the module information
     *dtype* is the kernel precision
+    *dim* is '1d' or '2d'
+    *result* is a vector to contain the results of the call
     The resulting call method takes the *pars*, a list of values for
 …
     Call :meth:`release` when done with the kernel instance.
     """
     def __init__(self, kernel, dtype, model_info, q_vectors):
+    def __init__(self, model, q_vectors):
         # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
+        q_input = GpuInput(q_vectors, dtype)
+        self.kernel = kernel
+        self.info = model_info
+        self.dtype = dtype
+        self.dim = '2d' if q_input.is_2d else '1d'
+        # plus three for the normalization values
+        self.result = np.empty(q_input.nq+1, dtype)
+        # Inputs and outputs for each kernel call
+        # Note: res may be shorter than res_b if global_size != nq
+        dtype = model.dtype
+        self.q_input = GpuInput(q_vectors, dtype)
+        self._model = model
+        # F16 isn't sufficient, so don't support it
+        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32
+        self._cache_key = unique_id()
+        # attributes accessed from the outside
+        self.dim = '2d' if self.q_input.is_2d else '1d'
+        self.info = model.info
+        self.dtype = model.dtype
+        # holding place for the returned value
+        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
+        extra_q = 4  # total weight, form volume, shell volume and R_eff
+        self.result = np.empty(self.q_input.nq*nout+extra_q, dtype)
+    @property
+    def _result_b(self):
+        """Lazy creation of result buffer so it can survive context reset"""
         env = environment()
+        self.queue = env.get_queue(dtype)
+        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
+                                  q_input.global_size[0] * dtype.itemsize)
+        self.q_input = q_input # allocated by GpuInput above
+        self._need_release = [self.result_b, self.q_input]
+        self.real = (np.float32 if dtype == generate.F32
+                     else np.float64 if dtype == generate.F64
+                     else np.float16 if dtype == generate.F16
+                     else np.float32)  # will never get here, so use np.float32
+    def __call__(self, call_details, values, cutoff, magnetic):
+        key = self._cache_key
+        if key not in env.cache:
+            context = env.context[self.dtype]
+            width = ((self.result.size+31)//32)*32 * self.dtype.itemsize
+            buffer = cl.Buffer(context, mf.READ_WRITE, width)
+            env.cache[key] = buffer
+        return env.cache[key]
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        context = self.queue.context
+        # Arrange data transfer to card
+        env = environment()
+        queue = env.queue[self._model.dtype]
+        context = queue.context
+        # Arrange data transfer to/from card
+        q_b = self.q_input.q_b
+        result_b = self._result_b
         details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
                               hostbuf=call_details.buffer)
 …
                              hostbuf=values)
+        kernel = self.kernel[1 if magnetic else 0]
+        args = [
+        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy'
+        kernel = self._model.fetch_kernel(name)
+        kernel_args = [
             np.uint32(self.q_input.nq), None, None,
+            details_b, values_b, self.q_input.q_b, self.result_b,
+            self.real(cutoff),
+            details_b, values_b, q_b, result_b,
+            self._as_dtype(cutoff),
+            np.uint32(effective_radius_type),
+        ]
         #print("Calling OpenCL")
         #call_details.show(values)
         # Call kernel and retrieve results
+        #Call kernel and retrieve results
         wait_for = None
         last_nap = time.clock()
+        last_nap = clock()
         step = 1000000//self.q_input.nq + 1
         for start in range(0, call_details.num_eval, step):
             stop = min(start + step, call_details.num_eval)
             #print("queuing",start,stop)
             args[1:3] = [np.int32(start), np.int32(stop)]
             wait_for = [kernel(self.queue, self.q_input.global_size, None,
                                *args, wait_for=wait_for)]
+            kernel_args[1:3] = [np.int32(start), np.int32(stop)]
+            wait_for = [kernel(queue, self.q_input.global_size, None,
+                               *kernel_args, wait_for=wait_for)]
             if stop < call_details.num_eval:
                 # Allow other processes to run
                 wait_for[0].wait()
                 current_time = time.clock()
+                current_time = clock()
                 if current_time - last_nap > 0.5:
                     time.sleep(0.05)
+                    time.sleep(0.001)
                     last_nap = current_time
         cl.enqueue_copy(self.queue, self.result, self.result_b)
+        cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for)
         #print("result", self.result)
 …
                 v.release()
-        pd_norm = self.result[self.q_input.nq]
-        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
-        background = values[1]
-        #print("scale",scale,values[0],self.result[self.q_input.nq],background)
-        return scale*self.result[:self.q_input.nq] + background
-        # return self.result[:self.q_input.nq]
     def release(self):
         # type: () -> None
 …
         Release resources associated with the kernel.
         """
+        for v in self._need_release:
+            v.release()
+        self._need_release = []
+        environment().free_buffer(id(self))
+        self.q_input.release()
     def __del__(self):

sasmodels/kerneldll.py

-                      r1a3559f
+                      re44432d
     pass
 # pylint: enable=unused-import
+# Compiler output is a byte stream that needs to be decode in python 3
+decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8'))
 if "SAS_DLL_PATH" in os.environ:
 …
         subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT)
     except subprocess.CalledProcessError as exc:
+        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output))
+        output = decode(exc.output)
+        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, output))
     if not os.path.exists(output):
         raise RuntimeError("compile failed.  File is in %r"%source)
 …
         # int, int, int, int*, double*, double*, double*, double*, double
         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type]
+        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32]
         names = [generate.kernel_name(self.info, variant)
                  for variant in ("Iq", "Iqxy", "Imagnetic")]
 …
     def __init__(self, kernel, model_info, q_input):
         # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None
+        #,model_info,q_input)
         self.kernel = kernel
         self.info = model_info
 …
         self.dtype = q_input.dtype
         self.dim = '2d' if q_input.is_2d else '1d'
+        self.result = np.empty(q_input.nq+1, q_input.dtype)
+        # leave room for f1/f2 results in case we need to compute beta for 1d models
+        nout = 2 if self.info.have_Fq else 1
+        # +4 for total weight, shell volume, effective radius, form volume
+        self.result = np.empty(q_input.nq*nout + 4, self.dtype)
         self.real = (np.float32 if self.q_input.dtype == generate.F32
                      else np.float64 if self.q_input.dtype == generate.F64
                      else np.float128)
+    def __call__(self, call_details, values, cutoff, magnetic):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
         kernel = self.kernel[1 if magnetic else 0]
         args = [
 …
             None, # pd_stop pd_stride[MAX_PD]
             call_details.buffer.ctypes.data, # problem
             values.ctypes.data,  #pars
             self.q_input.q.ctypes.data, #q
+            values.ctypes.data,  # pars
+            self.q_input.q.ctypes.data, # q
             self.result.ctypes.data,   # results
             self.real(cutoff), # cutoff
+            effective_radius_type, # cutoff
+        ]
         #print("Calling DLL")
 …
             kernel(*args) # type: ignore
-        #print("returned",self.q_input.q, self.result)
-        pd_norm = self.result[self.q_input.nq]
-        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
-        background = values[1]
-        #print("scale",scale,background)
-        return scale*self.result[:self.q_input.nq] + background
     def release(self):
         # type: () -> None

sasmodels/kernelpy.py

-                      r12eec1e
+                      raa8c6e0
 import numpy as np  # type: ignore
+from numpy import pi
+try:
+    from numpy import cbrt
+except ImportError:
+    def cbrt(x): return x ** (1.0/3.0)
 from .generate import F64
 …
         # Generate a closure which calls the form_volume if it exists.
+        form_volume = model_info.form_volume
+        self._volume = ((lambda: form_volume(*volume_args)) if form_volume else
+                        (lambda: 1.0))
+    def __call__(self, call_details, values, cutoff, magnetic):
+        self._volume_args = volume_args
+        volume = model_info.form_volume
+        shell = model_info.shell_volume
+        radius = model_info.effective_radius
+        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume
+                        else (lambda: [volume(*volume_args)]*2) if volume
+                        else (lambda: (1.0, 1.0)))
+        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius
+                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume
+                        else (lambda mode: 1.0))
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
         if magnetic:
 …
         #print("Calling python kernel")
         #call_details.show(values)
+        res = _loops(self._parameter_vector, self._form, self._volume,
+                     self.q_input.nq, call_details, values, cutoff)
+        return res
+        radius = ((lambda: 0.0) if effective_radius_type == 0
+                  else (lambda: self._radius(effective_radius_type)))
+        self.result = _loops(
+            self._parameter_vector, self._form, self._volume, radius,
+            self.q_input.nq, call_details, values, cutoff)
     def release(self):
 …
            form,          # type: Callable[[], np.ndarray]
            form_volume,   # type: Callable[[], float]
+           form_radius,   # type: Callable[[], float]
            nq,            # type: int
            call_details,  # type: details.CallDetails
            values,        # type: np.ndarray
            cutoff         # type: float
+           cutoff,        # type: float
           ):
     # type: (...) -> None
 …
     #                                                              #
     ################################################################
+    # WARNING: Trickery ahead
+    # The parameters[] vector is embedded in the closures for form(),
+    # form_volume() and form_radius().  We set the initial vector from
+    # the values for the model parameters. As we loop through the polydispesity
+    # mesh, we update the components with the polydispersity values before
+    # calling the respective functions.
     n_pars = len(parameters)
     parameters[:] = values[2:n_pars+2]
     if call_details.num_active == 0:
+        pd_norm = float(form_volume())
+        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
+        background = values[1]
+        return scale*form() + background
+    pd_value = values[2+n_pars:2+n_pars + call_details.num_weights]
+    pd_weight = values[2+n_pars + call_details.num_weights:]
+    pd_norm = 0.0
+    partial_weight = np.NaN
+    weight = np.NaN
+    p0_par = call_details.pd_par[0]
+    p0_length = call_details.pd_length[0]
+    p0_index = p0_length
+    p0_offset = call_details.pd_offset[0]
+    pd_par = call_details.pd_par[:call_details.num_active]
+    pd_offset = call_details.pd_offset[:call_details.num_active]
+    pd_stride = call_details.pd_stride[:call_details.num_active]
+    pd_length = call_details.pd_length[:call_details.num_active]
+    total = np.zeros(nq, 'd')
+    for loop_index in range(call_details.num_eval):
+        # update polydispersity parameter values
+        if p0_index == p0_length:
+            pd_index = (loop_index//pd_stride)%pd_length
+            parameters[pd_par] = pd_value[pd_offset+pd_index]
+            partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
+            p0_index = loop_index%p0_length
+        weight = partial_weight * pd_weight[p0_offset + p0_index]
+        parameters[p0_par] = pd_value[p0_offset + p0_index]
+        p0_index += 1
+        if weight > cutoff:
+            # Call the scattering function
+            # Assume that NaNs are only generated if the parameters are bad;
+            # exclude all q for that NaN.  Even better would be to have an
+            # INVALID expression like the C models, but that is too expensive.
+            Iq = np.asarray(form(), 'd')
+            if np.isnan(Iq).any():
+                continue
+            # update value and norm
+            total += weight * Iq
+            pd_norm += weight * form_volume()
+    scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
+    background = values[1]
+    return scale*total + background
+        total = form()
+        weight_norm = 1.0
+        weighted_shell, weighted_form = form_volume()
+        weighted_radius = form_radius()
+    else:
+        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights]
+        pd_weight = values[2+n_pars + call_details.num_weights:]
+        weight_norm = 0.0
+        weighted_form = 0.0
+        weighted_shell = 0.0
+        weighted_radius = 0.0
+        partial_weight = np.NaN
+        weight = np.NaN
+        p0_par = call_details.pd_par[0]
+        p0_length = call_details.pd_length[0]
+        p0_index = p0_length
+        p0_offset = call_details.pd_offset[0]
+        pd_par = call_details.pd_par[:call_details.num_active]
+        pd_offset = call_details.pd_offset[:call_details.num_active]
+        pd_stride = call_details.pd_stride[:call_details.num_active]
+        pd_length = call_details.pd_length[:call_details.num_active]
+        total = np.zeros(nq, 'd')
+        for loop_index in range(call_details.num_eval):
+            # update polydispersity parameter values
+            if p0_index == p0_length:
+                pd_index = (loop_index//pd_stride)%pd_length
+                parameters[pd_par] = pd_value[pd_offset+pd_index]
+                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
+                p0_index = loop_index%p0_length
+            weight = partial_weight * pd_weight[p0_offset + p0_index]
+            parameters[p0_par] = pd_value[p0_offset + p0_index]
+            p0_index += 1
+            if weight > cutoff:
+                # Call the scattering function
+                # Assume that NaNs are only generated if the parameters are bad;
+                # exclude all q for that NaN.  Even better would be to have an
+                # INVALID expression like the C models, but that is expensive.
+                Iq = np.asarray(form(), 'd')
+                if np.isnan(Iq).any():
+                    continue
+                # update value and norm
+                total += weight * Iq
+                weight_norm += weight
+                shell, form = form_volume()
+                weighted_form += weight * form
+                weighted_shell += weight * shell
+                weighted_radius += weight * form_radius()
+    result = np.hstack((total, weight_norm, weighted_form, weighted_shell, weighted_radius))
+    return result

sasmodels/mixture.py

recf895e	r39a06c9
143	143	#model_info.tests = []
144	144	#model_info.source = []
145		~~# Iq, Iqxy, form_volume, ER, VR and sesans~~
146	145	# Remember the component info blocks so we can build the model
147	146	model_info.composition = ('mixture', parts)

sasmodels/model_test.py

-                      r12eec1e
+                      r5024a56
 Usage::
     python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ...
+    python -m sasmodels.model_test [opencl|cuda|dll] model1 model2 ...
     if model1 is 'all', then all except the remaining models will be tested
 Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1),
 and the ER and VR are computed.  The return values at these points are not
+considered.  The test is only to verify that the models run to completion,
 and do not produce inf or NaN.
+and Fq is called to make sure R_eff, volume and volume ratio are computed.
+The return values at these points are not considered.  The test is only to
+verify that the models run to completion, and do not produce inf or NaN.
 Tests are defined with the *tests* attribute in the model.py file.  *tests*
 is a list of individual tests to run, where each test consists of the
 parameter values for the test, the q-values and the expected results.  For
+the effective radius test, the q-value should be 'ER'.  For the VR test,
+the q-value should be 'VR'.  For 1-D tests, either specify the q value or
+a list of q-values, and the corresponding I(q) value, or list of I(q) values.
+the effective radius test and volume ratio tests, use the extended output
+form, which checks each output of kernel.Fq. For 1-D tests, either specify
+the q value or a list of q-values, and the corresponding I(q) value, or
+list of I(q) values.
 That is::
 …
                         [I(qx1, qy1), I(qx2, qy2), ...]],
+        [ {parameters}, 'ER', ER(pars) ],
+        [ {parameters}, 'VR', VR(pars) ],
+        [ {parameters}, q, F(q), F^2(q), R_eff, V, V_r ],
         ...
+    ]
 …
 from . import core
 from .core import list_models, load_model_info, build_model
 from .direct_model import call_kernel, call_ER, call_VR
+from .direct_model import call_kernel, call_Fq
 from .exception import annotate_exception
 from .modelinfo import expand_pars
 from .kernelcl import use_opencl
+from .kernelcuda import use_cuda
+from . import product
 # pylint: disable=unused-import
 …
     Construct the pyunit test suite.
+    *loaders* is the list of kernel drivers to use, which is one of
+    *["dll", "opencl"]*, *["dll"]* or *["opencl"]*.  For python models,
+    the python driver is always used.
+    *loaders* is the list of kernel drivers to use (dll, opencl or cuda).
+    For python model the python driver is always used.
     *models* is the list of models to test, or *["all"]* to test all models.
 …
                                     test_method_name,
                                     platform="ocl", dtype=None,
+                                    stash=stash)
+            #print("defining", test_name)
+            suite.addTest(test)
+        # test using cuda if desired and available
+        if 'cuda' in loaders and use_cuda():
+            test_name = "%s-cuda"%model_name
+            test_method_name = "test_%s_cuda" % model_info.id
+            # Using dtype=None so that the models that are only
+            # correct for double precision are not tested using
+            # single precision.  The choice is determined by the
+            # presence of *single=False* in the model file.
+            test = ModelTestCase(test_name, model_info,
+                                    test_method_name,
+                                    platform="cuda", dtype=None,
                                     stash=stash)
             #print("defining", test_name)
 …
                 ({}, [0.001, 0.01, 0.1], [None]*3),
                 ({}, [(0.1, 0.1)]*2, [None]*2),
+                # test that ER/VR will run if they exist
+                ({}, 'ER', None),
+                ({}, 'VR', None),
+                # test that Fq will run, and return R_eff, V, V_r
+                ({}, 0.1, None, None, None, None, None),
+                ]
             tests = smoke_tests
 …
             if self.info.tests is not None:
                 tests += self.info.tests
+            S_tests = [test for test in tests if '@S' in test[0]]
+            P_tests = [test for test in tests if '@S' not in test[0]]
             try:
                 model = build_model(self.info, dtype=self.dtype,
                                     platform=self.platform)
+                results = [self.run_one(model, test) for test in tests]
+                results = [self.run_one(model, test) for test in P_tests]
+                for test in S_tests:
+                    # pull the S model name out of the test defn
+                    pars = test[0].copy()
+                    s_name = pars.pop('@S')
+                    ps_test = [pars] + list(test[1:])
+                    # build the P@S model
+                    s_info = load_model_info(s_name)
+                    ps_info = product.make_product_info(self.info, s_info)
+                    ps_model = build_model(ps_info, dtype=self.dtype,
+                                           platform=self.platform)
+                    # run the tests
+                    results.append(self.run_one(ps_model, ps_test))
                 if self.stash:
                     for test, target, actual in zip(tests, self.stash[0], results):
 …
                 # Check for missing tests.  Only do so for the "dll" tests
                 # to reduce noise from both opencl and dll, and because
+                # to reduce noise from both opencl and cuda, and because
                 # python kernels use platform="dll".
                 if self.platform == "dll":
 …
         def _find_missing_tests(self):
             # type: () -> None
+            """make sure there are 1D, 2D, ER and VR tests as appropriate"""
+            model_has_VR = callable(self.info.VR)
+            model_has_ER = callable(self.info.ER)
+            """make sure there are 1D and 2D tests as appropriate"""
             model_has_1D = True
             model_has_2D = any(p.type == 'orientation'
 …
             single = [test for test in self.info.tests
                       if not isinstance(test[2], list) and test[2] is not None]
-            tests_has_VR = any(test[1] == 'VR' for test in single)
-            tests_has_ER = any(test[1] == 'ER' for test in single)
             tests_has_1D_single = any(isinstance(test[1], float) for test in single)
             tests_has_2D_single = any(isinstance(test[1], tuple) for test in single)
 …
             missing = []
-            if model_has_VR and not tests_has_VR:
-                missing.append("VR")
-            if model_has_ER and not tests_has_ER:
-                missing.append("ER")
             if model_has_1D and not (tests_has_1D_single or tests_has_1D_multiple):
                 missing.append("1D")
 …
             # type: (KernelModel, TestCondition) -> None
             """Run a single test case."""
             user_pars, x, y = test
+            user_pars, x, y = test[:3]
             pars = expand_pars(self.info.parameters, user_pars)
             invalid = invalid_pars(self.info.parameters, pars)
 …
             self.assertEqual(len(y), len(x))
+            if x[0] == 'ER':
+                actual = np.array([call_ER(model.info, pars)])
+            elif x[0] == 'VR':
+                actual = np.array([call_VR(model.info, pars)])
+            elif isinstance(x[0], tuple):
+            if isinstance(x[0], tuple):
                 qx, qy = zip(*x)
                 q_vectors = [np.array(qx), np.array(qy)]
-                kernel = model.make_kernel(q_vectors)
-                actual = call_kernel(kernel, pars)
             else:
                 q_vectors = [np.array(x)]
+                kernel = model.make_kernel(q_vectors)
+            kernel = model.make_kernel(q_vectors)
+            if len(test) == 3:
                 actual = call_kernel(kernel, pars)
+            self.assertTrue(len(actual) > 0)
+            self.assertEqual(len(y), len(actual))
+            for xi, yi, actual_yi in zip(x, y, actual):
+                self._check_vectors(x, y, actual, 'I')
+                return actual
+            else:
+                y1 = y
+                y2 = test[3] if not isinstance(test[3], list) else [test[3]]
+                F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars)
+                if F1 is not None:  # F1 is none for models with Iq instead of Fq
+                    self._check_vectors(x, y1, F1, 'F')
+                self._check_vectors(x, y2, F2, 'F^2')
+                self._check_scalar(test[4], R_eff, 'R_eff')
+                self._check_scalar(test[5], volume, 'volume')
+                self._check_scalar(test[6], volume_ratio, 'form:shell ratio')
+                return F2
+        def _check_scalar(self, target, actual, name):
+            if target is None:
+                # smoke test --- make sure it runs and produces a value
+                self.assertTrue(not np.isnan(actual),
+                                'invalid %s: %s' % (name, actual))
+            elif np.isnan(target):
+                # make sure nans match
+                self.assertTrue(np.isnan(actual),
+                                '%s: expected:%s; actual:%s'
+                                % (name, target, actual))
+            else:
+                # is_near does not work for infinite values, so also test
+                # for exact values.
+                self.assertTrue(target == actual or is_near(target, actual, 5),
+                                '%s: expected:%s; actual:%s'
+                                % (name, target, actual))
+        def _check_vectors(self, x, target, actual, name='I'):
+            self.assertTrue(len(actual) > 0,
+                            '%s(...) expected return'%name)
+            if target is None:
+                return
+            self.assertEqual(len(target), len(actual),
+                             '%s(...) returned wrong length'%name)
+            for xi, yi, actual_yi in zip(x, target, actual):
                 if yi is None:
                     # smoke test --- make sure it runs and produces a value
                     self.assertTrue(not np.isnan(actual_yi),
                                     'invalid f(%s): %s' % (xi, actual_yi))
+                                    'invalid %s(%s): %s' % (name, xi, actual_yi))
                 elif np.isnan(yi):
+                    # make sure nans match
                     self.assertTrue(np.isnan(actual_yi),
                                     'f(%s): expected:%s; actual:%s'
                                     % (xi, yi, actual_yi))
+                                    '%s(%s): expected:%s; actual:%s'
+                                    % (name, xi, yi, actual_yi))
                 else:
                     # is_near does not work for infinite values, so also test
                     # for exact values.  Note that this will not
+                    # for exact values.
                     self.assertTrue(yi == actual_yi or is_near(yi, actual_yi, 5),
+                                    'f(%s); expected:%s; actual:%s'
+                                    % (xi, yi, actual_yi))
+            return actual
+                                    '%s(%s); expected:%s; actual:%s'
+                                    % (name, xi, yi, actual_yi))
     return ModelTestCase
 …
     invalid = []
     for par in sorted(pars.keys()):
+        # special handling of R_eff mode, which is not a usual parameter
+        if par == product.RADIUS_MODE_ID:
+            continue
         parts = par.split('_pd')
         if len(parts) > 1 and parts[1] not in ("", "_n", "nsigma", "type"):
 …
     # Build a test suite containing just the model
     loaders = ['opencl'] if use_opencl() else ['dll']
+    loaders = ['opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll']
     suite = unittest.TestSuite()
     _add_model_to_suite(loaders, suite, model_info)
 …
         loaders = ['opencl']
         models = models[1:]
+    elif models and models[0] == 'cuda':
+        if not use_cuda():
+            print("cuda is not available")
+            return 1
+        loaders = ['cuda']
+        models = models[1:]
     elif models and models[0] == 'dll':
         # TODO: test if compiler is available?
         loaders = ['dll']
         models = models[1:]
-    elif models and models[0] == 'opencl_and_dll':
-        loaders = ['opencl', 'dll'] if use_opencl() else ['dll']
-        models = models[1:]
     else:
+        loaders = ['opencl', 'dll'] if use_opencl() else ['dll']
+        loaders = ['dll']
+        if use_opencl():
+            loaders.append('opencl')
+        if use_cuda():
+            loaders.append('cuda')
     if not models:
         print("""\
 usage:
   python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ...
+  python -m sasmodels.model_test [-v] [opencl|cuda|dll] model1 model2 ...
 If -v is included on the command line, then use verbose output.
 If neither opencl nor dll is specified, then models will be tested with
 both OpenCL and dll; the compute target is ignored for pure python models.
+If no platform is specified, then models will be tested with dll, and
+if available, OpenCL and CUDA; the compute target is ignored for pure python models.
 If model1 is 'all', then all except the remaining models will be tested.
 …
     Run "nosetests sasmodels" on the command line to invoke it.
     """
+    loaders = ['opencl', 'dll'] if use_opencl() else ['dll']
+    loaders = ['dll']
+    if use_opencl():
+        loaders.append('opencl')
+    if use_cuda():
+        loaders.append('cuda')
     tests = make_suite(loaders, ['all'])
     def build_test(test):

sasmodels/modelinfo.py

-                      rbd547d0
+                      rc1799d3
     parameter.length = length
     parameter.length_control = control
     return parameter
 …
     will have a magnitude and a direction, which may be different from
     other sld parameters. The volume parameters are used for calls
     to form_volume within the kernel (required for volume normalization)
     and for calls to ER and VR for effective radius and volume ratio
     respectively.
+    to form_volume within the kernel (required for volume normalization),
+    to shell_volume (for hollow shapes), and to effective_radius (for
+    structure factor interactions) respectively.
     *description* is a short description of the parameter.  This will
 …
       parameters counted as n individual parameters p1, p2, ...
+    * *common_parameters* is the list of common parameters, with a unique
+      copy for each model so that structure factors can have a default
+      background of 0.0.
     * *call_parameters* is the complete list of parameters to the kernel,
       including scale and background, with vector parameters recorded as
 …
     parameters don't use vector notation, and instead use p1, p2, ...
     """
-    # scale and background are implicit parameters
-    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS]
     def __init__(self, parameters):
         # type: (List[Parameter]) -> None
+        # scale and background are implicit parameters
+        # Need them to be unique to each model in case they have different
+        # properties, such as default=0.0 for structure factor backgrounds.
+        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS]
         self.kernel_parameters = parameters
         self._set_vector_lengths()
         self.npars = sum(p.length for p in self.kernel_parameters)
         self.nmagnetic = sum(p.length for p in self.kernel_parameters
 …
         if self.nmagnetic:
             self.nvalues += 3 + 3*self.nmagnetic
         self.call_parameters = self._get_call_parameters()
         self.defaults = self._get_defaults()
 …
                          if p.polydisperse and p.type not in ('orientation', 'magnetic'))
         self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse)
+    def set_zero_background(self):
+        """
+        Set the default background to zero for this model.  This is done for
+        structure factor models.
+        """
+        # type: () -> None
+        # Make sure background is the second common parameter.
+        assert self.common_parameters[1].id == "background"
+        self.common_parameters[1].default = 0.0
+        self.defaults = self._get_defaults()
     def check_angles(self):
 …
     def _get_call_parameters(self):
         # type: () -> List[Parameter]
         full_list = self.COMMON[:]
+        full_list = self.common_parameters[:]
         for p in self.kernel_parameters:
             if p.length == 1:
 …
         # Gather the user parameters in order
         result = control + self.COMMON
+        result = control + self.common_parameters
         for p in self.kernel_parameters:
             if not is2d and p.type in ('orientation', 'magnetic'):
 …
 #: Set of variables defined in the model that might contain C code
 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code']
+C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code']
 def _find_source_lines(model_info, kernel_module):
 …
         # Custom sum/multi models
         return kernel_module.model_info
     info = ModelInfo()
+    # Build the parameter table
     #print("make parameter table", kernel_module.parameters)
     parameters = make_parameter_table(getattr(kernel_module, 'parameters', []))
+    # background defaults to zero for structure factor models
+    structure_factor = getattr(kernel_module, 'structure_factor', False)
+    if structure_factor:
+        parameters.set_zero_background()
+    # TODO: remove demo parameters
+    # The plots in the docs are generated from the model default values.
+    # Sascomp set parameters from the command line, and so doesn't need
+    # demo values for testing.
     demo = expand_pars(parameters, getattr(kernel_module, 'demo', None))
     filename = abspath(kernel_module.__file__).replace('.pyc', '.py')
     kernel_id = splitext(basename(filename))[0]
 …
     info.category = getattr(kernel_module, 'category', None)
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
+    # TODO: find Fq by inspection
+    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None)
+    info.have_Fq = getattr(kernel_module, 'have_Fq', False)
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
     # Note: custom.load_custom_kernel_module assumes the C sources are defined
 …
     info.source = getattr(kernel_module, 'source', [])
     info.c_code = getattr(kernel_module, 'c_code', None)
+    info.effective_radius = getattr(kernel_module, 'effective_radius', None)
     # TODO: check the structure of the tests
     info.tests = getattr(kernel_module, 'tests', [])
-    info.ER = getattr(kernel_module, 'ER', None) # type: ignore
-    info.VR = getattr(kernel_module, 'VR', None) # type: ignore
     info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore
+    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore
     info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore
     info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore
 …
     info.lineno = {}
     _find_source_lines(info, kernel_module)
     return info
 …
     #: provided in the file.
     structure_factor = None # type: bool
+    #: True if the model defines an Fq function with signature
+    #: void Fq(double q, double *F1, double *F2, ...)
+    have_Fq = False
+    #: List of options for computing the effective radius of the shape,
+    #: or None if the model is not usable as a form factor model.
+    effective_radius_type = None   # type: List[str]
     #: List of C source files used to define the model.  The source files
     #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the
     #: model defines orientation parameters. Files containing the most basic
     #: functions must appear first in the list, followed by the files that
+    #: use those functions.  Form factors are indicated by providing
+    #: an :attr:`ER` function.
+    #: use those functions.
     source = None           # type: List[str]
+    #: The set of tests that must pass.  The format of the tests is described
+    #: in :mod:`model_test`.
+    tests = None            # type: List[TestCondition]
+    #: Returns the effective radius of the model given its volume parameters.
+    #: The presence of *ER* indicates that the model is a form factor model
+    #: that may be used together with a structure factor to form an implicit
+    #: multiplication model.
+    #:
+    #: The parameters to the *ER* function must be marked with type *volume*.
+    #: in the parameter table.  They will appear in the same order as they
+    #: do in the table.  The values passed to *ER* will be vectors, with one
+    #: value for each polydispersity condition.  For example, if the model
+    #: is polydisperse over both length and radius, then both length and
+    #: radius will have the same number of values in the vector, with one
+    #: value for each *length X radius*.  If only *radius* is polydisperse,
+    #: then the value for *length* will be repeated once for each value of
+    #: *radius*.  The *ER* function should return one effective radius for
+    #: each parameter set.  Multiplicity parameters will be received as
+    #: arrays, with one row per polydispersity condition.
+    ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]]
+    #: Returns the occupied volume and the total volume for each parameter set.
+    #: See :attr:`ER` for details on the parameters.
+    VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]]
+    #: Arbitrary C code containing supporting functions, etc., to be inserted
+    #: after everything in source.  This can include Iq and Iqxy functions with
+    #: the full function signature, including all parameters.
+    c_code = None
+    #: inline source code, added after all elements of source
+    c_code = None           # type: Optional[str]
     #: Returns the form volume for python-based models.  Form volume is needed
     #: for volume normalization in the polydispersity integral.  If no
 …
     #: C code, either defined as a string, or in the sources.
     form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]]
+    #: Returns the shell volume for python-based models.  Form volume and
+    #: shell volume are needed for volume normalization in the polydispersity
+    #: integral and structure interactions for hollow shapes.  If no
+    #: parameters are *volume* parameters, then shell volume is not needed.
+    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq`
+    #: defined using a string containing C code), shell_volume must also be
+    #: C code, either defined as a string, or in the sources.
+    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]]
     #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined
     #: by the parameter table.  *Iq* can be defined as a python function, or
 …
     #: Line numbers for symbols defining C code
     lineno = None           # type: Dict[str, int]
+    #: The set of tests that must pass.  The format of the tests is described
+    #: in :mod:`model_test`.
+    tests = None            # type: List[TestCondition]
     def __init__(self):

sasmodels/models/_spherepy.py

-                      rca4444f
+                      r304c775
     return 1.333333333333333 * pi * radius ** 3
+def effective_radius(mode, radius):
+    return radius
 def Iq(q, sld, sld_solvent, radius):
     #print "q",q
 …
 sesans.vectorized = True  # sesans accepts an array of z values
-def ER(radius):
-    return radius
-# VR defaults to 1.0
 demo = dict(scale=1, background=0,
             sld=6, sld_solvent=1,

sasmodels/models/barbell.c

-                      r108e70e
+                      r99658f6
 static double
+Iq(double q, double sld, double solvent_sld,
+radius_from_excluded_volume(double radius_bell, double radius, double length)
+{
+    const double hdist = sqrt(square(radius_bell) - square(radius));
+    const double length_tot = length + 2.0*(hdist+ radius);
+    return 0.5*cbrt(0.75*radius_bell*(2.0*radius_bell*length_tot + (radius_bell + length_tot)*(M_PI*radius_bell + length_tot)));
+}
+static double
+radius_from_volume(double radius_bell, double radius, double length)
+{
+    const double vol_barbell = form_volume(radius_bell,radius,length);
+    return cbrt(vol_barbell/M_4PI_3);
+}
+static double
+radius_from_totallength(double radius_bell, double radius, double length)
+{
+    const double hdist = sqrt(square(radius_bell) - square(radius));
+    return 0.5*length + hdist + radius_bell;
+}
+static double
+effective_radius(int mode, double radius_bell, double radius, double length)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(radius_bell, radius , length);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(radius_bell, radius , length);
+    case 3: // radius
+        return radius;
+    case 4: // half length
+        return 0.5*length;
+    case 5: // half total length
+        return radius_from_totallength(radius_bell,radius,length);
+    }
+}
+static void
+Fq(double q,double *F1, double *F2, double sld, double solvent_sld,
     double radius_bell, double radius, double length)
+{
 …
     const double zm = M_PI_4;
     const double zb = M_PI_4;
+    double total = 0.0;
+    double total_F1 = 0.0;
+    double total_F2 = 0.0;
     for (int i = 0; i < GAUSS_N; i++){
         const double alpha= GAUSS_Z[i]*zm + zb;
 …
         SINCOS(alpha, sin_alpha, cos_alpha);
         const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length);
+        total += GAUSS_W[i] * Aq * Aq * sin_alpha;
+        total_F1 += GAUSS_W[i] * Aq * sin_alpha;
+        total_F2 += GAUSS_W[i] * Aq * Aq * sin_alpha;
+    }
     // translate dx in [-1,1] to dx in [lower,upper]
+    const double form = total*zm;
+    const double form_avg = total_F1*zm;
+    const double form_squared_avg = total_F2*zm;
     //Contrast
     const double s = (sld - solvent_sld);
+    return 1.0e-4 * s * s * form;
+    *F1 = 1.0e-2 * s * form_avg;
+    *F2 = 1.0e-4 * s * s * form_squared_avg;
+}
 static double

sasmodels/models/barbell.py

-                      r2d81cfe
+                      r99658f6
 .. [#] H Kaya and N R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda
    and errata)
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"]
+have_Fq = True
+effective_radius_type = [
+    "equivalent cylinder excluded volume","equivalent volume sphere", "radius", "half length", "half total length",
+    ]
 def random():

sasmodels/models/binary_hard_sphere.c

-                      r925ad6e
+                      r71b751d
 double form_volume(void);
 double Iq(double q,
+double Iq(double q,
     double lg_radius, double sm_radius,
     double lg_vol_frac, double sm_vol_frac,
     double lg_sld, double sm_sld, double solvent_sld
     );
 void calculate_psfs(double qval,
     double r2, double nf2,
 …
     double lg_vol_frac, double sm_vol_frac,
     double lg_sld, double sm_sld, double solvent_sld)
+{
     double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten;       //my local names
 …
     double phi1,phi2,phr,a3;
     double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
     r2 = lg_radius;
     r1 = sm_radius;
 …
     rho1 = sm_sld;
     rhos = solvent_sld;
     phi = phi1 + phi2;
     aa = r1/r2;
 …
     // calculate the PSF's here
     calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
     // /* do form factor calculations  */
     v1 = M_4PI_3*r1*r1*r1;
     v2 = M_4PI_3*r2*r2*r2;
     n1 = phi1/v1;
     n2 = phi2/v2;
     qr1 = r1*q;
     qr2 = r2*q;
 …
     inten *= 1.0e8;
     ///*convert rho^2 in 10^-6A to A*/
     inten *= 1.0e-12;
+    inten *= 1.0e-12;
     return(inten);
+}
 …
     double aa, double phi,
     double *s11, double *s22, double *s12)
+{
     //  variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
     //   calculate constant terms
     double s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
 …
     double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
     double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
     s2 = 2.0*r2;
 //    s1 = aa*s2;  why is this never used?  check original paper?
 …
     gm1=(v1*a1+a3*v2*a2)*.5;
     gm12=2.*gm1*(1.-aa)/aa;
     //c
+    //c
     //c   calculate the direct correlation functions and print results
     //c
     //  do 20 j=1,npts
     yy=qval*s2;
     //c   calculate direct correlation functions
 …
     t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
     f11=24.*v1*(t1+t2+t3)/ay3;
     //c ------c22
     y2=yy*yy;
 …
     tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
     f22=24.*v2*(tt1+tt2+tt3)/y3;
     //c   -----c12
     yl=.5*yy*(1.-aa);
 …
     ttt4=a1*(t41+t42)/y1;
     f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
     c11=f11;
     c22=f22;
     c12=f12;
     *s11=1./(1.+c11-(c12)*c12/(1.+c22));
     *s22=1./(1.+c22-(c12)*c12/(1.+c11));
     *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));
+    *s22=1./(1.+c22-(c12)*c12/(1.+c11));
+    *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));
     return;
+}

sasmodels/models/capped_cylinder.c

-                      r108e70e
+                      r99658f6
 static double
+Iq(double q, double sld, double solvent_sld,
+radius_from_excluded_volume(double radius, double radius_cap, double length)
+{
+    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius);
+    const double length_tot = length + 2.0*hc;
+    return 0.5*cbrt(0.75*radius*(2.0*radius*length_tot + (radius + length_tot)*(M_PI*radius + length_tot)));
+}
+static double
+radius_from_volume(double radius, double radius_cap, double length)
+{
+    const double vol_cappedcyl = form_volume(radius,radius_cap,length);
+    return cbrt(vol_cappedcyl/M_4PI_3);
+}
+static double
+radius_from_totallength(double radius, double radius_cap, double length)
+{
+    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius);
+    return 0.5*length + hc;
+}
+static double
+effective_radius(int mode, double radius, double radius_cap, double length)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(radius, radius_cap, length);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(radius, radius_cap, length);
+    case 3: // radius
+        return radius;
+    case 4: // half length
+        return 0.5*length;
+    case 5: // half total length
+        return radius_from_totallength(radius, radius_cap,length);
+    }
+}
+static void
+Fq(double q,double *F1, double *F2, double sld, double solvent_sld,
     double radius, double radius_cap, double length)
+{
 …
     const double zm = M_PI_4;
     const double zb = M_PI_4;
+    double total = 0.0;
+    double total_F1 = 0.0;
+    double total_F2 = 0.0;
     for (int i=0; i<GAUSS_N ;i++) {
         const double theta = GAUSS_Z[i]*zm + zb;
 …
         const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
         // scale by sin_theta for spherical coord integration
+        total += GAUSS_W[i] * Aq * Aq * sin_theta;
+        total_F1 += GAUSS_W[i] * Aq * sin_theta;
+        total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta;
+    }
     // translate dx in [-1,1] to dx in [lower,upper]
+    const double form = total * zm;
+    const double form_avg = total_F1 * zm;
+    const double form_squared_avg = total_F2 * zm;
     // Contrast
     const double s = (sld - solvent_sld);
+    return 1.0e-4 * s * s * form;
+    *F1 = 1.0e-2 * s * form_avg;
+    *F2 = 1.0e-4 * s * s * form_squared_avg;
+}

sasmodels/models/capped_cylinder.py

-                      r2d81cfe
+                      r99658f6
 .. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda
    and errata)
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"]
+have_Fq = True
+effective_radius_type = [
+    "equivalent cylinder excluded volume", "equivalent volume sphere", "radius", "half length", "half total length",
+    ]
 def random():

sasmodels/models/core_multi_shell.c

-                      rc3ccaec
+                      rd42dd4a
 static double
 form_volume(double core_radius, double fp_n, double thickness[])
+outer_radius(double core_radius, double fp_n, double thickness[])
+{
   double r = core_radius;
 …
     r += thickness[i];
+  }
   return M_4PI_3 * cube(r);
+  return r;
+}
 static double
+Iq(double q, double core_sld, double core_radius,
+form_volume(double core_radius, double fp_n, double thickness[])
+{
+  return M_4PI_3 * cube(outer_radius(core_radius, fp_n, thickness));
+}
+static double
+effective_radius(int mode, double core_radius, double fp_n, double thickness[])
+{
+  switch (mode) {
+  default:
+  case 1: // outer radius
+    return outer_radius(core_radius, fp_n, thickness);
+  case 2: // core radius
+    return core_radius;
+  }
+}
+static void
+Fq(double q, double *F1, double *F2, double core_sld, double core_radius,
    double solvent_sld, double fp_n, double sld[], double thickness[])
+{
 …
+  }
   f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sas_3j1x_x(q*r);
+  return f * f * 1.0e-4;
+  *F1 = 1e-2 * f;
+  *F2 = 1e-4 * f * f;
+}

sasmodels/models/core_multi_shell.py

-                      r2d81cfe
+                      ree60aa7
 source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"]
+have_Fq = True
+effective_radius_type = ["outer radius", "core radius"]
 def random():
 …
     return np.asarray(z), np.asarray(rho)
-def ER(radius, n, thickness):
-    """Effective radius"""
-    n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse
-    return np.sum(thickness[:n], axis=0) + radius
 demo = dict(sld_core=6.4,
             radius=60,

sasmodels/models/core_shell_bicelle.c

-                      r108e70e
+                      r99658f6
 static double
+Iq(double q,
+radius_from_excluded_volume(double radius, double thick_rim, double thick_face, double length)
+{
+    const double radius_tot = radius + thick_rim;
+    const double length_tot = length + 2.0*thick_face;
+    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot)));
+}
+static double
+radius_from_volume(double radius, double thick_rim, double thick_face, double length)
+{
+    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length);
+    return cbrt(volume_bicelle/M_4PI_3);
+}
+static double
+radius_from_diagonal(double radius, double thick_rim, double thick_face, double length)
+{
+    const double radius_tot = radius + thick_rim;
+    const double length_tot = length + 2.0*thick_face;
+    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot);
+}
+static double
+effective_radius(int mode, double radius, double thick_rim, double thick_face, double length)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(radius, thick_rim, thick_face, length);
+    case 2: // equivalent sphere
+        return radius_from_volume(radius, thick_rim, thick_face, length);
+    case 3: // outer rim radius
+        return radius + thick_rim;
+    case 4: // half outer thickness
+        return 0.5*length + thick_face;
+    case 5: // half diagonal
+        return radius_from_diagonal(radius,thick_rim,thick_face,length);
+    }
+}
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
     double radius,
     double thick_radius,
 …
     const double halflength = 0.5*length;
+    double total = 0.0;
+    double total_F1 = 0.0;
+    double total_F2 = 0.0;
     for(int i=0;i<GAUSS_N;i++) {
         double theta = (GAUSS_Z[i] + 1.0)*uplim;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
         double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
+        double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
                                    halflength, sld_core, sld_face, sld_rim, sld_solvent);
+        total += GAUSS_W[i]*fq*fq*sin_theta;
+        total_F1 += GAUSS_W[i]*form*sin_theta;
+        total_F2 += GAUSS_W[i]*form*form*sin_theta;
+    }
+    // Correct for integration range
+    total_F1 *= uplim;
+    total_F2 *= uplim;
+    // calculate value of integral to return
+    double answer = total*uplim;
+    return 1.0e-4*answer;
+    *F1 = 1.0e-2*total_F1;
+    *F2 = 1.0e-4*total_F2;
+}

sasmodels/models/core_shell_bicelle.py

-                      r2d81cfe
+                      r99658f6
    from Proquest <http://search.proquest.com/docview/304915826?accountid
    =26379>`_
+   L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",
           "core_shell_bicelle.c"]
+have_Fq = True
+effective_radius_type = [
+    "excluded volume","equivalent volume sphere", "outer rim radius",
+    "half outer thickness", "half diagonal",
+    ]
 def random():

sasmodels/models/core_shell_bicelle_elliptical.c

-                      r108e70e
+                      r99658f6
 static double
+Iq(double q,
+radius_from_excluded_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double r_equiv     = sqrt((r_minor + thick_rim)*(r_minor*x_core + thick_rim));
+    const double length_tot  = length + 2.0*thick_face;
+    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_tot + (r_equiv + length_tot)*(M_PI*r_equiv + length_tot)));
+}
+static double
+radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length);
+    return cbrt(volume_bicelle/M_4PI_3);
+}
+static double
+radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor);
+    const double radius_max_tot = radius_max + thick_rim;
+    const double length_tot = length + 2.0*thick_face;
+    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot);
+}
+static double
+effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(r_minor, x_core, thick_rim, thick_face, length);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);
+    case 3: // outer rim average radius
+        return 0.5*r_minor*(1.0 + x_core) + thick_rim;
+    case 4: // outer rim min radius
+        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 5: // outer max radius
+        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 6: // half outer thickness
+        return 0.5*length + thick_face;
+    case 7: // half diagonal
+        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length);
+    }
+}
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
     double r_minor,
     double x_core,
 …
     //initialize integral
+    double outer_sum = 0.0;
+    double outer_total_F1 = 0.0;
+    double outer_total_F2 = 0.0;
     for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
-        //const double va = 0.0;
-        //const double vb = 1.0;
         //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
         const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0;
 …
         const double si1 = sas_sinx_x(halfheight*qc);
         const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
+        double inner_sum=0.0;
+        double inner_total_F1 = 0;
+        double inner_total_F2 = 0;
         for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
+            // inner integral limits
+            //const double vaj=0.0;
+            //const double vbj=M_PI;
+            //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2;
+            const double rr = sqrt(r2A - r2B*cos(phi));
+            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
+            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2;
+            const double rr = sqrt(r2A - r2B*cos(beta));
             const double be1 = sas_2J1x_x(rr*qab);
             const double be2 = sas_2J1x_x((rr+thick_rim)*qab);
             const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
+            const double f = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1;
+            inner_sum += GAUSS_W[j] * fq * fq;
+            inner_total_F1 += GAUSS_W[j] * f;
+            inner_total_F2 += GAUSS_W[j] * f * f;
+        }
         //now calculate outer integral
+        outer_sum += GAUSS_W[i] * inner_sum;
+        outer_total_F1 += GAUSS_W[i] * inner_total_F1;
+        outer_total_F2 += GAUSS_W[i] * inner_total_F2;
+    }
+    // now complete change of integration variables (1-0)/(1-(-1))= 0.5
+    outer_total_F1 *= 0.25;
+    outer_total_F2 *= 0.25;
+    return outer_sum*2.5e-05;
+    //convert from [1e-12 A-1] to [cm-1]
+    *F1 = 1e-2*outer_total_F1;
+    *F2 = 1e-4*outer_total_F2;
+}

sasmodels/models/core_shell_bicelle_elliptical.py

-                      rfc3ae1b
+                      r99658f6
 .. [#]
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",
           "core_shell_bicelle_elliptical.c"]
+have_Fq = True
+effective_radius_type = [
+    "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius",
+    "outer max radius", "half outer thickness", "half diagonal",
+    ]
 def random():
 …
 tests = [
     [{'radius': 30.0, 'x_core': 3.0,
       'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'ER', 1],
     [{'radius': 30.0, 'x_core': 3.0,
       'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'VR', 1],
+    #[{'radius': 30.0, 'x_core': 3.0,
+    #  'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'ER', 1],
+    #[{'radius': 30.0, 'x_core': 3.0,
+    #  'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'VR', 1],
     [{'radius': 30.0, 'x_core': 3.0,

sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

-                      r108e70e
+                      r99658f6
 static double
+Iq(double q,
+radius_from_excluded_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double r_equiv     = sqrt((r_minor + thick_rim)*(r_minor*x_core + thick_rim));
+    const double length_tot  = length + 2.0*thick_face;
+    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_tot + (r_equiv + length_tot)*(M_PI*r_equiv + length_tot)));
+}
+static double
+radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length);
+    return cbrt(volume_bicelle/M_4PI_3);
+}
+static double
+radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor);
+    const double radius_max_tot = radius_max + thick_rim;
+    const double length_tot = length + 2.0*thick_face;
+    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot);
+}
+static double
+effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(r_minor, x_core, thick_rim, thick_face, length);
+    case 2: // equivalent sphere
+        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);
+    case 3: // outer rim average radius
+        return 0.5*r_minor*(1.0 + x_core) + thick_rim;
+    case 4: // outer rim min radius
+        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 5: // outer max radius
+        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);
+    case 6: // half outer thickness
+        return 0.5*length + thick_face;
+    case 7: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face
+            //  or thick_rim is extremely large)
+        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length);
+    }
+}
+static void
+Fq(double q,
+        double *F1,
+        double *F2,
         double r_minor,
         double x_core,
 …
         double sigma)
+{
-    double si1,si2,be1,be2;
      // core_shell_bicelle_elliptical_belt, RKH 5th Oct 2017, core_shell_bicelle_elliptical
      // tested briefly against limiting cases of cylinder, hollow cylinder & elliptical cylinder models
 …
     const double r2A = 0.5*(square(r_major) + square(r_minor));
     const double r2B = 0.5*(square(r_major) - square(r_minor));
+    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
+    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*halfheight;
+    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
     // dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face,
+    const double dr1 = (-rhor - rhoh + rhoc + rhosolv) *M_PI*r_minor*r_major*
+.0*halfheight;
+    const double dr2 = (rhor-rhosolv) *M_PI*(r_minor+thick_rim)*(
+         r_major+thick_rim)* 2.0*halfheight;
+    const double dr3 = (rhoh-rhosolv) *M_PI*r_minor*r_major*
+.0*(halfheight+thick_face);
+    const double dr1 = vol1*(-rhor - rhoh + rhoc + rhosolv);
+    const double dr2 = vol2*(rhor-rhosolv);
+    const double dr3 = vol3*(rhoh-rhosolv);
     //initialize integral
+    double outer_sum = 0.0;
+    double outer_total_F1 = 0.0;
+    double outer_total_F2 = 0.0;
     for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         // since we generate these lots of times, why not store them somewhere?
+        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0;
+        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
+        double inner_sum=0;
+        double sinarg1 = q*halfheight*cos_alpha;
+        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
+        si1 = sas_sinx_x(sinarg1);
+        si2 = sas_sinx_x(sinarg2);
+        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0;
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double si1 = sas_sinx_x(halfheight*qc);
+        const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
+        double inner_total_F1 = 0;
+        double inner_total_F2 = 0;
         for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
 …
             const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(beta));
+            double besarg1 = q*rr*sin_alpha;
+            double besarg2 = q*(rr+thick_rim)*sin_alpha;
+            be1 = sas_2J1x_x(besarg1);
+            be2 = sas_2J1x_x(besarg2);
+            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 +
+                                              dr2*si1*be2 +
+                                              dr3*si2*be1);
+            const double be1 = sas_2J1x_x(rr*qab);
+            const double be2 = sas_2J1x_x((rr+thick_rim)*qab);
+            const double f = dr1*si1*be1 + dr2*si1*be2 + dr3*si2*be1;
+            inner_total_F1 += GAUSS_W[j] * f;
+            inner_total_F2 += GAUSS_W[j] * f * f;
+        }
         //now calculate outer integral
+        outer_sum += GAUSS_W[i] * inner_sum;
+        outer_total_F1 += GAUSS_W[i] * inner_total_F1;
+        outer_total_F2 += GAUSS_W[i] * inner_total_F2;
+    }
+    // now complete change of integration variables (1-0)/(1-(-1))= 0.5
+    outer_total_F1 *= 0.25;
+    outer_total_F2 *= 0.25;
+    return outer_sum*2.5e-05*exp(-0.5*square(q*sigma));
+    //convert from [1e-12 A-1] to [cm-1]
+    *F1 = 1e-2*outer_total_F1*exp(-0.25*square(q*sigma));
+    *F2 = 1e-4*outer_total_F2*exp(-0.5*square(q*sigma));
+}
 …
     const double si1 = sas_sinx_x( halfheight*qc );
     const double si2 = sas_sinx_x( (halfheight + thick_face)*qc );
+    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1);
+    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma));
+    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1;
+    const double atten = exp(-0.5*(qa*qa + qb*qb + qc*qc)*(sigma*sigma));
+    return 1.0e-4 * fq*fq * atten;
+}

sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

-                      rfc3ae1b
+                      r99658f6
 .. [#]
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",
           "core_shell_bicelle_elliptical_belt_rough.c"]
+have_Fq = True
+effective_radius_type = [
+    "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius",
+    "outer max radius", "half outer thickness", "half diagonal",
+    ]
 demo = dict(scale=1, background=0,
 …
 tests = [
     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'ER', 1],
     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1],
+    #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'ER', 1],
+    #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1],
     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0,

sasmodels/models/core_shell_cylinder.c

-                      rd86f0fc
+                      r99658f6
 static double
+Iq(double q,
+radius_from_excluded_volume(double radius, double thickness, double length)
+{
+    const double radius_tot = radius + thickness;
+    const double length_tot = length + 2.0*thickness;
+    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot)));
+}
+static double
+radius_from_volume(double radius, double thickness, double length)
+{
+    const double volume_outer_cyl = form_volume(radius,thickness,length);
+    return cbrt(volume_outer_cyl/M_4PI_3);
+}
+static double
+radius_from_diagonal(double radius, double thickness, double length)
+{
+    const double radius_outer = radius + thickness;
+    const double length_outer = length + 2.0*thickness;
+    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer);
+}
+static double
+effective_radius(int mode, double radius, double thickness, double length)
+{
+    switch (mode) {
+    default:
+    case 1: //cylinder excluded volume
+        return radius_from_excluded_volume(radius, thickness, length);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(radius, thickness, length);
+    case 3: // outer radius
+        return radius + thickness;
+    case 4: // half outer length
+        return 0.5*length + thickness;
+    case 5: // half min outer length
+        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness);
+    case 6: // half max outer length
+        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness);
+    case 7: // half outer diagonal
+        return radius_from_diagonal(radius,thickness,length);
+    }
+}
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
     double core_sld,
     double shell_sld,
 …
     const double shell_h = (0.5*length + thickness);
     const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
+    double total = 0.0;
+    double total_F1 = 0.0;
+    double total_F2 = 0.0;
     for (int i=0; i<GAUSS_N ;i++) {
         // translate a point in [-1,1] to a point in [0, pi/2]
 …
         const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
             + _cyl(shell_vd, shell_r*qab, shell_h*qc);
+        total += GAUSS_W[i] * fq * fq * sin_theta;
+        total_F1 += GAUSS_W[i] * fq * sin_theta;
+        total_F2 += GAUSS_W[i] * fq * fq * sin_theta;
+    }
     // translate dx in [-1,1] to dx in [lower,upper]
     //const double form = (upper-lower)/2.0*total;
+    return 1.0e-4 * total * M_PI_4;
+    *F1 = 1.0e-2 * total_F1 * M_PI_4;
+    *F2 = 1.0e-4 * total_F2 * M_PI_4;
+}
 static double

sasmodels/models/core_shell_cylinder.py

-                      re31b19a
+                      r99658f6
 -1452
 .. [#kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"]
+def ER(radius, thickness, length):
+    """
+    Returns the effective radius used in the S*P calculation
+    """
+    radius = radius + thickness
+    length = length + 2 * thickness
+    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
+    return 0.5 * (ddd) ** (1. / 3.)
+have_Fq = True
+effective_radius_type = [
+    "excluded volume", "equivalent volume sphere", "outer radius", "half outer length",
+    "half min outer dimension", "half max outer dimension", "half outer diagonal",
+    ]
 def random():

sasmodels/models/core_shell_ellipsoid.c

-                      r108e70e
+                      r99658f6
 static double
+Iq(double q,
+radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    return cbrt(volume_ellipsoid/M_4PI_3);
+}
+static double
+radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    // Trivial cases
+    if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell;
+    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.;
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double radius_equat_tot = radius_equat_core + thick_shell;
+    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
+    const double ratio = (radius_polar_tot < radius_equat_tot
+                          ? radius_polar_tot / radius_equat_tot
+                          : radius_equat_tot / radius_polar_tot);
+    const double e1 = sqrt(1.0 - ratio*ratio);
+    const double b1 = 1.0 + asin(e1) / (e1 * ratio);
+    const double bL = (1.0 + e1) / (1.0 - e1);
+    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL);
+    const double delta = 0.75 * b1 * b2;
+    const double ddd = 2.0 * (delta + 1.0) * radius_polar_tot * radius_equat_tot * radius_equat_tot;
+    return 0.5 * cbrt(ddd);
+}
+static double
+effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)
+{
+    const double radius_equat_tot = radius_equat_core + thick_shell;
+    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell;
+    switch (mode) {
+    default:
+    case 1: // average outer curvature
+        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);
+    case 3: // min outer radius
+        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot);
+    case 4: // max outer radius
+        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot);
+    }
+}
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
     double radius_equat_core,
     double x_core,
 …
     const double m = 0.5;
     const double b = 0.5;
+    double total = 0.0;     //initialize intergral
+    double total_F1 = 0.0;     //initialize intergral
+    double total_F2 = 0.0;     //initialize intergral
     for(int i=0;i<GAUSS_N;i++) {
         const double cos_theta = GAUSS_Z[i]*m + b;
 …
             equat_shell, polar_shell,
             sld_core_shell, sld_shell_solvent);
+        total += GAUSS_W[i] * fq * fq;
+        total_F1 += GAUSS_W[i] * fq;
+        total_F2 += GAUSS_W[i] * fq * fq;
+    }
+    total *= m;
+    total_F1 *= m;
+    total_F2 *= m;
     // convert to [cm-1]
+    return 1.0e-4 * total;
+    *F1 = 1.0e-2 * total_F1;
+    *F2 = 1.0e-4 * total_F2;
+}
 static double

sasmodels/models/core_shell_ellipsoid.py

-                      r2d81cfe
+                      r99658f6
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]
+def ER(radius_equat_core, x_core, thick_shell, x_polar_shell):
+    """
+        Returns the effective radius used in the S*P calculation
+    """
+    from .ellipsoid import ER as ellipsoid_ER
+    polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell
+    equat_outer = radius_equat_core + thick_shell
+    return ellipsoid_ER(polar_outer, equat_outer)
+have_Fq = True
+effective_radius_type = [
+    "average outer curvature", "equivalent volume sphere",
+    "min outer radius", "max outer radius",
+    ]
 def random():

sasmodels/models/core_shell_parallelepiped.c

-                      rdbf1a60
+                      r99658f6
 static double
+Iq(double q,
+radius_from_excluded_volume(double length_a, double length_b, double length_c,
+                   double thick_rim_a, double thick_rim_b, double thick_rim_c)
+{
+    double r_equiv, length;
+    double lengths[3] = {length_a+thick_rim_a, length_b+thick_rim_b, length_c+thick_rim_c};
+    double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2]));
+    double length_1 = lengthmax;
+    double length_2 = lengthmax;
+    double length_3 = lengthmax;
+    for(int ilen=0; ilen<3; ilen++) {
+        if (lengths[ilen] < length_1) {
+            length_2 = length_1;
+            length_1 = lengths[ilen];
+            } else {
+                if (lengths[ilen] < length_2) {
+                        length_2 = lengths[ilen];
+                }
+            }
+    }
+    if(length_2-length_1 > length_3-length_2) {
+        r_equiv = sqrt(length_2*length_3/M_PI);
+        length  = length_1;
+    } else  {
+        r_equiv = sqrt(length_1*length_2/M_PI);
+        length  = length_3;
+    }
+    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length)));
+}
+static double
+radius_from_volume(double length_a, double length_b, double length_c,
+                   double thick_rim_a, double thick_rim_b, double thick_rim_c)
+{
+    const double volume = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    return cbrt(volume/M_4PI_3);
+}
+static double
+radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b)
+{
+    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a;
+    return sqrt(area_xsec_paral/M_PI);
+}
+static double
+effective_radius(int mode, double length_a, double length_b, double length_c,
+                 double thick_rim_a, double thick_rim_b, double thick_rim_c)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    case 3: // half outer length a
+        return 0.5 * length_a + thick_rim_a;
+    case 4: // half outer length b
+        return 0.5 * length_b + thick_rim_b;
+    case 5: // half outer length c
+        return 0.5 * length_c + thick_rim_c;
+    case 6: // equivalent circular cross-section
+        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b);
+    case 7: // half outer ab diagonal
+        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b));
+    case 8: // half outer diagonal
+        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c));
+    }
+}
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
     double core_sld,
     double arim_sld,
 …
     // outer integral (with gauss points), integration limits = 0, 1
     // substitute d_cos_alpha for sin_alpha d_alpha
+    double outer_sum = 0; //initialize integral
+    double outer_sum_F1 = 0; //initialize integral
+    double outer_sum_F2 = 0; //initialize integral
     for( int i=0; i<GAUSS_N; i++) {
         const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
 …
         // inner integral (with gauss points), integration limits = 0, 1
         // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta)
+        double inner_sum = 0.0;
+        double inner_sum_F1 = 0.0;
+        double inner_sum_F2 = 0.0;
         for(int j=0; j<GAUSS_N; j++) {
             const double u = 0.5 * ( GAUSS_Z[j] + 1.0 );
 …
 #endif
+            inner_sum += GAUSS_W[j] * f * f;
+            inner_sum_F1 += GAUSS_W[j] * f;
+            inner_sum_F2 += GAUSS_W[j] * f * f;
+        }
         // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
         inner_sum *= 0.5;
         // now sum up the outer integral
         outer_sum += GAUSS_W[i] * inner_sum;
+        // and sum up the outer integral
+        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * 0.5;
+        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * 0.5;
+    }
     // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
+    outer_sum *= 0.5;
+    outer_sum_F1 *= 0.5;
+    outer_sum_F2 *= 0.5;
     //convert from [1e-12 A-1] to [cm-1]
+    return 1.0e-4 * outer_sum;
+    *F1 = 1.0e-2 * outer_sum_F1;
+    *F2 = 1.0e-4 * outer_sum_F2;
+}

sasmodels/models/core_shell_parallelepiped.py

-                      rf89ec96
+                      r99658f6
 is the scattering length of the solvent.
 .. note::
+.. note::
    the code actually implements two substitutions: $d(cos\alpha)$ is
 …
    from Proquest <http://search.proquest.com/docview/304915826?accountid
    =26379>`_
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/gauss76.c", "core_shell_parallelepiped.c"]
+def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c):
+    """
+        Return equivalent radius (ER)
+    """
+    from .parallelepiped import ER as ER_p
+    a = length_a + 2*thick_rim_a
+    b = length_b + 2*thick_rim_b
+    c = length_c + 2*thick_rim_c
+    return ER_p(a, b, c)
+# VR defaults to 1.0
+have_Fq = True
+effective_radius_type = [
+    "equivalent cylinder excluded volume",
+    "equivalent volume sphere",
+    "half outer length_a", "half outer length_b", "half outer length_c",
+    "equivalent circular cross-section",
+    "half outer ab diagonal", "half outer diagonal",
+    ]
 def random():

sasmodels/models/core_shell_sphere.c

-                      r3a48772
+                      rd42dd4a
+double form_volume(double radius, double thickness);
+double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld);
+static double
+form_volume(double radius, double thickness)
+{
+    return M_4PI_3 * cube(radius + thickness);
+}
+double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) {
+static double
+effective_radius(int mode, double radius, double thickness)
+{
+    switch (mode) {
+    default:
+    case 1: // outer radius
+        return radius + thickness;
+    case 2: // core radius
+        return radius;
+    }
+}
+    double intensity = core_shell_kernel(q,
+static void
+Fq(double q, double *F1, double *F2, double radius,
+   double thickness, double core_sld, double shell_sld, double solvent_sld) {
+    double form = core_shell_fq(q,
                               radius,
                               thickness,
 …
                               shell_sld,
                               solvent_sld);
+    return intensity;
+    *F1 = 1.0e-2*form;
+    *F2 = 1.0e-4*form*form;
+}
-double form_volume(double radius, double thickness)
+{
-    return M_4PI_3 * cube(radius + thickness);
+}

sasmodels/models/core_shell_sphere.py

-                      rda1c8d1
+                      r304c775
 title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure."
 description = """
     F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3
                    + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
+    F(q) = [V_c (sld_core-sld_shell) 3 (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3
+            + V_s (sld_shell-sld_solvent) 3 (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3]
             V_s: Volume of the sphere shell
 …
 source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"]
+have_Fq = True
+effective_radius_type = ["outer radius", "core radius"]
 demo = dict(scale=1, background=0, radius=60, thickness=10,
             sld_core=1.0, sld_shell=2.0, sld_solvent=0.0)
-def ER(radius, thickness):
-    """
-        Equivalent radius
-        @param radius: core radius
-        @param thickness: shell thickness
-    """
-    return radius + thickness
 def random():
 …
 tests = [
     [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
+    [{'radius': 20.0, 'thickness': 10.0}, 0.1, None, None, 30.0, 4.*pi/3*30**3, 1.0],
     # The SasView test result was 0.00169, with a background of 0.001

sasmodels/models/cylinder.c

-                      r108e70e
+                      r99658f6
 static double
 fq(double qab, double qc, double radius, double length)
+_fq(double qab, double qc, double radius, double length)
+{
     return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
 …
 static double
+orient_avg_1D(double q, double radius, double length)
+radius_from_excluded_volume(double radius, double length)
+{
+    return 0.5*cbrt(0.75*radius*(2.0*radius*length + (radius + length)*(M_PI*radius + length)));
+}
+static double
+radius_from_volume(double radius, double length)
+{
+    return cbrt(M_PI*radius*radius*length/M_4PI_3);
+}
+static double
+radius_from_diagonal(double radius, double length)
+{
+    return sqrt(radius*radius + 0.25*length*length);
+}
+static double
+effective_radius(int mode, double radius, double length)
+{
+    switch (mode) {
+    default:
+    case 1:
+        return radius_from_excluded_volume(radius, length);
+    case 2:
+        return radius_from_volume(radius, length);
+    case 3:
+        return radius;
+    case 4:
+        return 0.5*length;
+    case 5:
+        return (radius < 0.5*length ? radius : 0.5*length);
+    case 6:
+        return (radius > 0.5*length ? radius : 0.5*length);
+    case 7:
+        return radius_from_diagonal(radius,length);
+    }
+}
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
+    double sld,
+    double solvent_sld,
+    double radius,
+    double length)
+{
     // translate a point in [-1,1] to a point in [0, pi/2]
 …
     const double zb = M_PI_4;
+    double total = 0.0;
+    double total_F1 = 0.0;
+    double total_F2 = 0.0;
     for (int i=0; i<GAUSS_N ;i++) {
         const double theta = GAUSS_Z[i]*zm + zb;
 …
         // theta (theta,phi) the projection of the cylinder on the detector plane
         SINCOS(theta , sin_theta, cos_theta);
+        const double form = fq(q*sin_theta, q*cos_theta, radius, length);
+        total += GAUSS_W[i] * form * form * sin_theta;
+        const double form = _fq(q*sin_theta, q*cos_theta, radius, length);
+        total_F1 += GAUSS_W[i] * form * sin_theta;
+        total_F2 += GAUSS_W[i] * form * form * sin_theta;
+    }
     // translate dx in [-1,1] to dx in [lower,upper]
+    return total*zm;
+    total_F1 *= zm;
+    total_F2 *= zm;
+    const double s = (sld - solvent_sld) * form_volume(radius, length);
+    *F1 = 1e-2 * s * total_F1;
+    *F2 = 1e-4 * s * s * total_F2;
+}
+static double
+Iq(double q,
+    double sld,
+    double solvent_sld,
+    double radius,
+    double length)
+{
+    const double s = (sld - solvent_sld) * form_volume(radius, length);
+    return 1.0e-4 * s * s * orient_avg_1D(q, radius, length);
+}
 static double
 …
     double length)
+{
+    const double form = _fq(qab, qc, radius, length);
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    const double form = fq(qab, qc, radius, length);
     return 1.0e-4 * square(s * form);
+}

sasmodels/models/cylinder.py

-                      r2d81cfe
+                      r99658f6
 J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997).
 G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951).
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 """
 …
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"]
+def ER(radius, length):
+    """
+        Return equivalent radius (ER)
+    """
+    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius))
+    return 0.5 * (ddd) ** (1. / 3.)
+have_Fq = True
+effective_radius_type = [
+    "excluded volume", "equivalent volume sphere", "radius",
+    "half length", "half min dimension", "half max dimension", "half diagonal",
+    ]
 def random():
 …
             phi_pd=10, phi_pd_n=5)
+# pylint: disable=bad-whitespace, line-too-long
 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
 # After redefinition of angles, find new tests values.  Was 10 10 in old coords
 …
+]
 del qx, qy  # not necessary to delete, but cleaner
+# Default radius and length
+radius, length = parameters[2][2], parameters[3][2]
+tests.extend([
+    ({'radius_effective_mode': 0}, 0.1, None, None, 0., pi*radius**2*length, 1.0),
+    ({'radius_effective_mode': 1}, 0.1, None, None, 0.5*(0.75*radius*(2.0*radius*length + (radius + length)*(pi*radius + length)))**(1./3.), None, None),
+    ({'radius_effective_mode': 2}, 0.1, None, None, (0.75*radius**2*length)**(1./3.), None, None),
+    ({'radius_effective_mode': 3}, 0.1, None, None, radius, None, None),
+    ({'radius_effective_mode': 4}, 0.1, None, None, length/2., None, None),
+    ({'radius_effective_mode': 5}, 0.1, None, None, min(radius, length/2.), None, None),
+    ({'radius_effective_mode': 6}, 0.1, None, None, max(radius, length/2.), None, None),
+    ({'radius_effective_mode': 7}, 0.1, None, None, np.sqrt(4*radius**2 + length**2)/2., None, None),
+])
+del radius, length
+# pylint: enable=bad-whitespace, line-too-long
 # ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc

sasmodels/models/dab.py

-                      r2d81cfe
+                      r304c775
     return pars
-# ER defaults to 1.0
-# VR defaults to 1.0
 demo = dict(scale=1, background=0, cor_length=50)

sasmodels/models/ellipsoid.c

-                      r108e70e
+                      r99658f6
+}
+static  double
+Iq(double q,
+static double
+radius_from_volume(double radius_polar, double radius_equatorial)
+{
+    return cbrt(radius_polar*radius_equatorial*radius_equatorial);
+}
+static double
+radius_from_curvature(double radius_polar, double radius_equatorial)
+{
+    // Trivial cases
+    if (radius_polar == radius_equatorial) return radius_polar;
+    if (radius_polar * radius_equatorial == 0.)  return 0.;
+    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    const double ratio = (radius_polar < radius_equatorial
+                          ? radius_polar / radius_equatorial
+                          : radius_equatorial / radius_polar);
+    const double e1 = sqrt(1.0 - ratio*ratio);
+    const double b1 = 1.0 + asin(e1) / (e1 * ratio);
+    const double bL = (1.0 + e1) / (1.0 - e1);
+    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL);
+    const double delta = 0.75 * b1 * b2;
+    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial;
+    return 0.5 * cbrt(ddd);
+}
+static double
+effective_radius(int mode, double radius_polar, double radius_equatorial)
+{
+    switch (mode) {
+    default:
+    case 1: // average curvature
+        return radius_from_curvature(radius_polar, radius_equatorial);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(radius_polar, radius_equatorial);
+    case 3: // min radius
+        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial);
+    case 4: // max radius
+        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial);
+    }
+}
+static void
+Fq(double q,
+    double *F1,
+    double *F2,
     double sld,
     double sld_solvent,
 …
     //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du
     const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0;
     // translate a point in [-1,1] to a point in [0, 1]
     // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;
     const double zm = 0.5;
     const double zb = 0.5;
+    double total = 0.0;
+    double total_F2 = 0.0;
+    double total_F1 = 0.0;
     for (int i=0;i<GAUSS_N;i++) {
         const double u = GAUSS_Z[i]*zm + zb;
         const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one);
         const double f = sas_3j1x_x(q*r);
+        total += GAUSS_W[i] * f * f;
+        total_F2 += GAUSS_W[i] * f * f;
+        total_F1 += GAUSS_W[i] * f;
+    }
     // translate dx in [-1,1] to dx in [lower,upper]
+    const double form = total*zm;
+    total_F1 *= zm;
+    total_F2 *= zm;
     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
+    return 1.0e-4 * s * s * form;
+    *F1 = 1e-2 * s * total_F1;
+    *F2 = 1e-4 * s * s * total_F2;
+}

sasmodels/models/ellipsoid.py

-                      r0168844
+                      r99658f6
+             ]
 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"]
+def ER(radius_polar, radius_equatorial):
+    # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    ee = np.empty_like(radius_polar)
+    idx = radius_polar > radius_equatorial
+    ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2
+    idx = radius_polar < radius_equatorial
+    ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2
+    valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial)
+    bd = 1.0 - ee[valid]
+    e1 = np.sqrt(ee[valid])
+    b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd))
+    bL = (1.0 + e1) / (1.0 - e1)
+    b2 = 1.0 + bd / 2 / e1 * np.log(bL)
+    delta = 0.75 * b1 * b2
+    ddd = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid]
+    r = np.zeros_like(radius_polar)
+    r[valid] = 0.5 * cbrt(ddd)
+    idx = radius_polar == radius_equatorial
+    r[idx] = radius_polar[idx]
+    return r
+have_Fq = True
+effective_radius_type = [
+    "average curvature", "equivalent volume sphere", "min radius", "max radius",
+    ]
 def random():

sasmodels/models/elliptical_cylinder.c

-                      r108e70e
+                      r99658f6
 static double
+Iq(double q, double radius_minor, double r_ratio, double length,
+radius_from_excluded_volume(double radius_minor, double r_ratio, double length)
+{
+    const double r_equiv = sqrt(radius_minor*radius_minor*r_ratio);
+    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length)));
+}
+static double
+radius_from_volume(double radius_minor, double r_ratio, double length)
+{
+    const double volume_ellcyl = form_volume(radius_minor,r_ratio,length);
+    return cbrt(volume_ellcyl/M_4PI_3);
+}
+static double
+radius_from_min_dimension(double radius_minor, double r_ratio, double hlength)
+{
+    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
+    return (rad_min < hlength ? rad_min : hlength);
+}
+static double
+radius_from_max_dimension(double radius_minor, double r_ratio, double hlength)
+{
+    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
+    return (rad_max > hlength ? rad_max : hlength);
+}
+static double
+radius_from_diagonal(double radius_minor, double r_ratio, double length)
+{
+    const double radius_max = (r_ratio > 1.0 ? radius_minor*r_ratio : radius_minor);
+    return sqrt(radius_max*radius_max + 0.25*length*length);
+}
+static double
+effective_radius(int mode, double radius_minor, double r_ratio, double length)
+{
+    switch (mode) {
+    default:
+    case 1: // equivalent cylinder excluded volume
+        return radius_from_excluded_volume(radius_minor, r_ratio, length);
+    case 2: // equivalent volume sphere
+        return radius_from_volume(radius_minor, r_ratio, length);
+    case 3: // average radius
+        return 0.5*radius_minor*(1.0 + r_ratio);
+    case 4: // min radius
+        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
+    case 5: // max radius
+        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
+    case 6: // equivalent circular cross-section
+        return sqrt(radius_minor*radius_minor*r_ratio);
+    case 7: // half length
+        return 0.5*length;
+    case 8: // half min dimension
+        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length);
+    case 9: // half max dimension
+        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length);
+    case 10: // half diagonal
+        return radius_from_diagonal(radius_minor,r_ratio,length);
+    }
+}
+static void
+Fq(double q, double *F1, double *F2, double radius_minor, double r_ratio, double length,
    double sld, double solvent_sld)
+{
 …
     //initialize integral
+    double outer_sum = 0.0;
+    double outer_sum_F1 = 0.0;
+    double outer_sum_F2 = 0.0;
     for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
 …
         const double sin_val = sqrt(1.0 - cos_val*cos_val);
         //const double arg = radius_minor*sin_val;
+        double inner_sum=0;
+        double inner_sum_F1 = 0.0;
+        double inner_sum_F2 = 0.0;
         for(int j=0;j<GAUSS_N;j++) {
             const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             const double r = sin_val*sqrt(rA - rB*cos(theta));
             const double be = sas_2J1x_x(q*r);
+            inner_sum += GAUSS_W[j] * be * be;
+            inner_sum_F1 += GAUSS_W[j] * be;
+            inner_sum_F2 += GAUSS_W[j] * be * be;
+        }
         //now calculate the value of the inner integral
+        inner_sum *= 0.5*(vbj-vaj);
+        inner_sum_F1 *= 0.5*(vbj-vaj);
+        inner_sum_F2 *= 0.5*(vbj-vaj);
         //now calculate outer integral
         const double si = sas_sinx_x(q*0.5*length*cos_val);
+        outer_sum += GAUSS_W[i] * inner_sum * si * si;
+        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * si;
+        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * si * si;
+    }
+    outer_sum *= 0.5*(vb-va);
+    //divide integral by Pi
+    const double form = outer_sum/M_PI;
+    // correct limits and divide integral by pi
+    outer_sum_F1 *= 0.5*(vb-va)/M_PI;
+    outer_sum_F2 *= 0.5*(vb-va)/M_PI;
     // scale by contrast and volume, and convert to to 1/cm units
+    const double vol = form_volume(radius_minor, r_ratio, length);
+    const double delrho = sld - solvent_sld;
+    return 1.0e-4*square(delrho*vol)*form;
+    const double volume = form_volume(radius_minor, r_ratio, length);
+    const double contrast = sld - solvent_sld;
+    const double s = contrast*volume;
+    *F1 = 1.0e-2*s*outer_sum_F1;
+    *F2 = 1.0e-4*s*s*outer_sum_F2;
+}
 …
     const double be = sas_2J1x_x(qr);
     const double si = sas_sinx_x(qc*0.5*length);
     const double Aq = be * si;
     const double delrho = sld - solvent_sld;
     const double vol = form_volume(radius_minor, r_ratio, length);
     return 1.0e-4 * square(delrho * vol * Aq);
+    const double fq = be * si;
+    const double contrast = sld - solvent_sld;
+    const double volume = form_volume(radius_minor, r_ratio, length);
+    return 1.0e-4 * square(contrast * volume * fq);
+}

sasmodels/models/elliptical_cylinder.py

-                      ra261a83
+                      r99658f6
 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and
 Neutron Scattering*, Plenum, New York, (1987) [see table 3.4]
+L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).
 Authorship and Verification
 …
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"]
+have_Fq = True
+effective_radius_type = [
+    "equivalent cylinder excluded volume", "equivalent volume sphere", "average radius", "min radius", "max radius",
+    "equivalent circular cross-section",
+    "half length", "half min dimension", "half max dimension", "half diagonal",
+    ]
 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0,
             sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30,
             theta_pd=10, phi_pd=2, psi_pd=3)
-def ER(radius_minor, axis_ratio, length):
-    """
-        Equivalent radius
-        @param radius_minor: Ellipse minor radius
-        @param axis_ratio: Ratio of major radius over minor radius
-        @param length: Length of the cylinder
-    """
-    radius = sqrt(radius_minor * radius_minor * axis_ratio)
-    ddd = 0.75 * radius * (2 * radius * length
-                           + (length + radius) * (length + pi * radius))
-    return 0.5 * (ddd) ** (1. / 3.)
 def random():
 …
 tests = [
     [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],
     [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1],
+#    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],
+#    [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1],
     # The SasView test result was 0.00169, with a background of 0.001

sasmodels/models/fcc_paracrystal.c

r108e70e	r71b751d
79	79	}
80	80
81
82	81	static double Iqabc(double qa, double qb, double qc,
83	82	double dnn, double d_factor, double radius,

sasmodels/models/flexible_cylinder_elliptical.c

r74768cb	r71b751d
72	72	return result;
73	73	}
74

sasmodels/models/fractal.c

-                      r4788822
+                      r71b751d
     //calculate P(q) for the spherical subunits
     const double pq = square(form_volume(radius) * (sld_block-sld_solvent)
                       *sas_3j1x_x(q*radius));
+    const double fq = form_volume(radius) * (sld_block-sld_solvent)
+                      *sas_3j1x_x(q*radius);
     // scale to units cm-1 sr-1 (assuming data on absolute scale)
 …
     //    combined: 1e-4 * I(q)
     return 1.e-4 * volfraction * sq * pq;
+    return 1.e-4 * volfraction * sq * fq * fq;
+}

sasmodels/models/fractal_core_shell.c

-                      rbdd08df
+                      rbad3093
    double cor_length)
+{
     //The radius for the building block of the core shell particle that is
+    //The radius for the building block of the core shell particle that is
     //needed by the Teixeira fractal S(q) is the radius of the whole particle.
     const double cs_radius = radius + thickness;
     const double sq = fractal_sq(q, cs_radius, fractal_dim, cor_length);
     const double pq = core_shell_kernel(q, radius, thickness,
                                         core_sld, shell_sld, solvent_sld);
+    const double fq = core_shell_fq(q, radius, thickness,
+                                    core_sld, shell_sld, solvent_sld);
+    // Note: core_shell_kernel already performs the 1e-4 unit conversion
+    return volfraction * sq * pq;
+    return 1.0e-4 * volfraction * sq * fq * fq;
+}

sasmodels/models/fractal_core_shell.py

-                      reb3eb38
+                      r2cc8aa2
     return radius + thickness
 tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
+#tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
+tests = [
 #         # At some point the SasView 3.x test result was deemed incorrect. The
           #following tests were verified against NIST IGOR macros ver 7.850.

sasmodels/models/fuzzy_sphere.py

-                      r2d81cfe
+                      ree60aa7
               ["sld_solvent", "1e-6/Ang^2",  3, [-inf, inf], "sld",    "Solvent scattering length density"],
               ["radius",      "Ang",        60, [0, inf],    "volume", "Sphere radius"],
               ["fuzziness",   "Ang",        10, [0, inf],    "",       "std deviation of Gaussian convolution for interface (must be << radius)"],
+              ["fuzziness",   "Ang",        10, [0, inf],    "volume",       "std deviation of Gaussian convolution for interface (must be << radius)"],
+             ]
 # pylint: enable=bad-whitespace,line-too-long
+source = ["lib/sas_3j1x_x.c"]
+# No volume normalization despite having a volume parameter
+# This should perhaps be volume normalized?
+form_volume = """
+    return M_4PI_3*cube(radius);
+    """
+Iq = """
+    const double qr = q*radius;
+    const double bes = sas_3j1x_x(qr);
+    const double qf = q*fuzziness;
+    const double fq = bes * (sld - sld_solvent) * form_volume(radius) * exp(-0.5*qf*qf);
+    return 1.0e-4*fq*fq;
+    """
+def ER(radius):
+    """
+    Return radius
+    """
+    return radius
+# VR defaults to 1.0
+source = ["lib/sas_3j1x_x.c","fuzzy_sphere.c"]
+have_Fq = True
+effective_radius_type = ["radius", "radius + fuzziness"]
 def random():

sasmodels/models/gaussian_peak.py

r2d81cfe	r304c775
50	50	"""
51	51
52		~~# VR defaults to 1.0~~
53
54	52	def random():
55	53	peak_pos = 10**np.random.uniform(-3, -1)

sasmodels/models/hardsphere.py

-                      r2d81cfe
+                      rc1799d3
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["radius_effective", "Ang", 50.0, [0, inf], "volume",
+parameters = [["radius_effective", "Ang", 50.0, [0, inf], "",
                "effective radius of hard sphere"],
               ["volfraction", "", 0.2, [0, 0.74], "",
                "volume fraction of hard spheres"],
+             ]
-# No volume normalization despite having a volume parameter
-# This should perhaps be volume normalized?
-form_volume = """
-    return 1.0;
-    """
 Iq = r"""
 …
     return pars
-# ER defaults to 0.0
-# VR defaults to 1.0
-demo = dict(radius_effective=200, volfraction=0.2,
-            radius_effective_pd=0.1, radius_effective_pd_n=40)
 # Q=0.001 is in the Taylor series, low Q part, so add Q=0.1,
 # assuming double precision sasview is correct

sasmodels/models/hayter_msa.py

r2d81cfe	r304c775
89	89	return 1.0;
90	90	"""
91		~~# ER defaults to 0.0~~
92		~~# VR defaults to 1.0~~
93	91
94	92	def random():

sasmodels/models/hollow_cylinder.c

-                      r108e70e
+                      r99658f6
 //#define INVALID(v) (v.radius_core >= v.radius)
-// From Igor library
 static double
 _hollow_cylinder_scaling(double integrand, double delrho, double volume)
+shell_volume(double radius, double thickness, double length)
+{
+    return 1.0e-4 * square(volume * delrho) * integrand;
+    return M_PI*length*(square(radius+thickness) - radius*radius);
+}
+static double
+form_volume(double radius, double thickness, double length)
+{
+    return M_PI*length*square(radius+thickness);
+}
+static double
+radius_from_excluded_volume(double radius, double thickness, double length)
+{

Changeset 9150036 in sasmodels

Legend: