Changeset 41917b8 – SasView

.travis.yml

-                      rb419c2d
+                      rce8c388
 language: python
+sudo:  false
+sudo: false
 matrix:
   include:
+    - os: linux
+      env:
+        - PY=2.7
+    - os: linux
+      env:
+        - PY=3.6
+    - os: osx
+      language: generic
+      env:
+        - PY=2.7
+    - os: osx
+      language: generic
+      env:
+        - PY=3.5
+# whitelist
+  - os: linux
+    env:
+    - PY=2.7
+  - os: linux
+    env:
+    - PY=3.6
+  - os: osx
+    language: generic
+    env:
+    - PY=2.7
+  - os: osx
+    language: generic
+    env:
+    - PY=3.5
 branches:
   only:
+    - master
+  - master
 addons:
   apt:
+    packages:
+      opencl-headers
+    packages: opencl-headers
 before_install:
+  - echo $TRAVIS_OS_NAME
+  - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then
+      wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
+    fi;
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then
+      wget https://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O miniconda.sh;
+    fi;
+  - bash miniconda.sh -b -p $HOME/miniconda
+  - export PATH="$HOME/miniconda/bin:$PATH"
+  - hash -r
+  - conda update --yes conda
+  # Useful for debugging any issues with conda
+  - conda info -a
+  - conda install --yes python=$PY numpy scipy cython mako cffi
+  # Not testing with opencl below, so don't need to install it
+  #- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then
+  #    pip install pyopencl;
+  #  fi;
+- if [[ $encrypted_cb04388797b6_iv ]]; then openssl aes-256-cbc -K $encrypted_cb04388797b6_key -iv $encrypted_cb04388797b6_iv
+  -in .travis/travis_rsa.enc -out .travis/travis_rsa -d; fi;
+- echo $TRAVIS_OS_NAME
+- if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
+  -O miniconda.sh; fi;
+- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then wget https://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh
+  -O miniconda.sh; fi;
+- bash miniconda.sh -b -p $HOME/miniconda
+- export PATH="$HOME/miniconda/bin:$PATH"
+- hash -r
+- conda update --yes conda
+- conda info -a
+- conda install --yes python=$PY numpy scipy cython mako cffi
 install:
+  - pip install bumps
+  - pip install unittest-xml-reporting
+- pip install bumps
+- pip install unittest-xml-reporting
 script:
 - python --version
 - python -m sasmodels.model_test -v dll all
+before_deploy:
+- echo -e "Host danse.chem.utk.edu\n\tStrictHostKeyChecking no\n" >> ~/.ssh/config
+deploy:
+  skip_cleanup: true
+  provider: script
+  script: "/bin/sh -ex ./deploy.sh"
+  on:
+    branch: master
 notifications:
   slack:

doc/conf.py

r8ae8532	r30b60d2
211	211	#latex_preamble = ''
212	212	LATEX_PREAMBLE=r"""
	213	\newcommand{\lt}{<}
	214	\newcommand{\gt}{>}
213	215	\renewcommand{\AA}{\text{\r{A}}} % Allow \AA in math mode
214	216	\usepackage[utf8]{inputenc} % Allow unicode symbols in text

doc/genmodel.py

r745b7bb	rb866abf
3	3	import sys, os, math, re
4	4	import numpy as np
	5	import matplotlib
	6	matplotlib.use('Agg')
5	7	import matplotlib.pyplot as plt
6	8	sys.path.insert(0, os.path.abspath('..'))

doc/guide/index.rst

r2e66ef5	rc0d7ab3
9	9	intro.rst
10	10	install.rst
	11	gpu_setup.rst
11	12	pd/polydispersity.rst
12	13	resolution.rst

doc/guide/install.rst

-                      rf8a2baa
+                      rc0d7ab3
 This will allow you to edit the files in the package and have the changes
 show up immediately in python the next time you load your program.
-OpenCL Installation
-*******************
-*Warning! GPU devices do not in general offer the same level of memory
-protection as CPU devices. If your code attempts to write outside allocated
-memory buffers unpredicatable behaviour may result (eg, your video display
-may freeze, or your system may crash, etc). Do not install OpenCL drivers
-without first checking for known issues (eg, some computer manufacturers
-install modified graphics drivers so replacing these may not be a good
-idea!). If in doubt, seek advice from an IT professional before proceeding
-further.*
-Check if you have OpenCL already installed
-==========================================
-**Windows**
-The following instructions are based on
-http://web.engr.oregonstate.edu/~mjb/cs475/DoIHaveOpenCL.pdf
-* Go to: Start -> Control Panel -> System & Security -> Administrative Tools
-* Double Click on Computer Managment
-* Click on Device Manager
-* Click open Display Adapters
-* Right-click on available adapter and select Properties
-* Click on Driver
-* Go to Driver Details
-* Scroll down and see if OpenCL is installed (look for OpenCL*.dll files)
-**Mac OSX**
-For OS X operating systems higher than 10.6 OpenCL is shipped along with
-the system.
-However, OpenCL has had a rocky history on Macs. Apple provide a useful
-compatibility table at https://support.apple.com/en-us/HT202823
-Installation
-============
-**Windows**
-Depending on the graphic card in your system, drivers
-can be obtained from different sources:
-* NVIDIA: https://developer.nvidia.com/opencl
-* AMD: http://developer.amd.com/tools-and-sdks/opencl-zone/
-**Mac OSX**
-N/A
-You cannot download OpenCL driver updates for your Mac. They are packaged
-with the normal quarterly OS X updates from Apple.
-.. note::
-    Intel provides OpenCL drivers for Intel processors at
-    https://software.intel.com/en-us/articles/opencl-drivers
-    These can sometimes make use of special vector instructions across multiple
-    processors, so it is worth installing if the GPU does not support double
-    precision. You can install this driver alongside the GPU driver for NVIDIA
-    or AMD.
-GPU Selection
-*************
-sasmodels evaluations can run on your graphics card (GPU) or they can run
-on the processor (CPU). In general, calculations performed on the GPU will run faster.
-To run on the GPU, your computer must have OpenCL drivers installed.
-For information about OpenCL installation see this
-:ref:`opencl-installation` guidance.
-Where the model is evaluated is a little bit complicated.
-If the model has the line *single=False* then it requires double precision.
-If the GPU is single precision only, then it will try running via OpenCL
-on the CPU.  If the OpenCL driver is not available for the CPU then
-it will run as a normal program on the CPU.
-For models with a large number of parameters or with a lot of code,
-the GPU may be too small to run the program effectively.
-In this case, you should try simplifying the model, maybe breaking it
-into several different modules so that you don't need *IF* statements in your code.
-If it is still too big, you can set *opencl=False* in the model file and
-the model will only run as a normal program on the CPU.
-This will not usually be necessary.
-Device Selection
-================
-If you have multiple GPU devices you can tell SasView which device to use.
-By default, SasView looks for one GPU and one CPU device
-from available OpenCL platforms.
-SasView prefers AMD or NVIDIA drivers for GPU, and prefers Intel or
-Apple drivers for CPU. Both GPU and CPU are included on the assumption that CPU
-is always available and supports double precision.
-The device order is important: GPU is checked before CPU on the assumption that
-it will be faster. By examining ~/sasview.log you can see which device SasView
-chose 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
-from the python console in SasView::
-    import pyopencl as cl
-    cl.create_some_context()
-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*.
-Compiler Selection
-==================
-For models run as normal programs, you may need to specify a compiler.
-This is done using the SAS_COMPILER environment variable.
-Set it to *tinycc* for the tinycc compiler, *msvc* for the
-Microsoft Visual C compiler, or *mingw* for the MinGW compiler.
-TinyCC is provided with SasView so that is the default.
-If you want one of the other compilers, be sure to have it available
-in your *PATH* so SasView can find it!
-*Document History*
-| 2017-05-17 Paul Kienzle

doc/guide/magnetism/magnetism.rst

-                      r990d8df
+                      r1f058ea
 .. figure::
     mag_img/mag_vector.bmp
+    mag_img/mag_vector.png
 The magnetic scattering length density is then
 …
 .. figure::
     mag_img/M_angles_pic.bmp
+    mag_img/M_angles_pic.png
 If the angles of the $Q$ vector and the spin-axis $x'$ to the $x$ - axis are
 …
 ===========   ================================================================
  M0_sld        = $D_M M_0$
  Up_theta      = $\theta_{up}$
+ Up_theta      = $\theta_\mathrm{up}$
  M_theta       = $\theta_M$
  M_phi         = $\phi_M$
 …
     The values of the 'Up_frac_i' and 'Up_frac_f' must be in the range 0 to 1.
+.. note::
+    This help document was last changed by Steve King, 02May2015
+*Document History*
+* Document History *
+| 2015-05-02 Steve King
 | 2017-05-08 Paul Kienzle

doc/guide/pd/polydispersity.rst

-                      rf8a2baa
+                      r1f058ea
            \exp\left(-\frac{(x - \bar x)^2}{2\sigma^2}\right)
 where $\bar x$ is the mean of the distribution and *Norm* is a normalization factor
 which is determined during the numerical calculation.
+where $\bar x$ is the mean of the distribution and *Norm* is a normalization
+factor which is determined during the numerical calculation.
 The polydispersity is
 …
 during the numerical calculation.
 The median value for the distribution will be the value given for the respective
 size parameter, for example, *radius=60*.
+The median value for the distribution will be the value given for the
+respective size parameter, for example, *radius=60*.
 The polydispersity is given by $\sigma$
 …
 Many commercial Dynamic Light Scattering (DLS) instruments produce a size
 polydispersity parameter, sometimes even given the symbol $p$ This
+polydispersity parameter, sometimes even given the symbol $p$\ ! This
 parameter is defined as the relative standard deviation coefficient of
 variation of the size distribution and is NOT the same as the polydispersity

doc/guide/plugin.rst

-                      r870a2f4
+                      r3048ec6
      at the start of the model documentation and as a tooltip in the SasView GUI
 - a **short discription**:
+- a **short description**:
    - this is the **description** string in the *.py* file
    - this is a medium length description which appears when you click
 …
 Models that do not conform to these requirements will *never* be incorporated
 into the built-in library.
-More complete documentation for the sasmodels package can be found at
-`<http://www.sasview.org/sasmodels>`_. In particular,
-`<http://www.sasview.org/sasmodels/api/generate.html#module-sasmodels.generate>`_
-describes the structure of a model.
 …
 **name = "mymodel"** defines the name of the model that is shown to the user.
+If it is not provided, it will use the name of the model file, with '_'
+replaced by spaces and the parts capitalized.  So *adsorbed_layer.py* will
+become *Adsorbed Layer*.  The predefined models all use the name of the
+model file as the name of the model, so the default may be changed.
+If it is not provided it will use the name of the model file. The name must
+be a valid variable name, starting with a letter and contains only letters,
+numbers or underscore.  Spaces, dashes, and other symbols are not permitted.
 **title = "short description"** is short description of the model which
 …
 - **"name"** is the name of the parameter shown on the FitPage.
+  - the name must be a valid variable name, starting with a letter and
+    containing only letters, numbers and underscore.
   - parameter names should follow the mathematical convention; e.g.,
     *radius_core* not *core_radius*, or *sld_solvent* not *solvent_sld*.
 …
     dispersion.
+Some models will have integer parameters, such as number of pearls in the
+pearl necklace model, or number of shells in the multi-layer vesicle model.
+The optimizers in BUMPS treat all parameters as floating point numbers which
+can take arbitrary values, even for integer parameters, so your model should
+round the incoming parameter value to the nearest integer inside your model
+you should round to the nearest integer.  In C code, you can do this using::
+    static double
+    Iq(double q, ..., double fp_n, ...)
+    {
+        int n = (int)(fp_n + 0.5);
+        ...
+    }
+in python::
+    def Iq(q, ..., n, ...):
+        n = int(n+0.5)
+        ...
+Derivative based optimizers such as Levenberg-Marquardt will not work
+for integer parameters since the partial derivative is always zero, but
+the remaining optimizers (DREAM, differential evolution, Nelder-Mead simplex)
+will still function.
 Model Computation
 …
      \approx \frac{\sum_{i=1}^n G(x_i) I(q,x_i)}{\sum_{i=1}^n G(x_i) V(x_i)}
 That is, the indivdual models do not need to include polydispersity
+That is, the individual models do not need to include polydispersity
 calculations, but instead rely on numerical integration to compute the
 appropriately smeared pattern.   Angular dispersion values over polar angle
 …
 The parameters *par1, par2, ...* are the list of non-orientation parameters
 to the model in the order that they appear in the parameter table.
 **Note that the autogenerated model file uses** *x* **rather than** *q*.
+**Note that the auto-generated model file uses** *x* **rather than** *q*.
 The *.py* file should import trigonometric and exponential functions from
 …
     sas_gamma(x):
         Gamma function $\text{sas_gamma}(x) = \Gamma(x)$.
+        Gamma function sas_gamma\ $(x) = \Gamma(x)$.
         The standard math function, tgamma(x) is unstable for $x < 1$
 …
     sas_erf(x), sas_erfc(x):
         Error function
         $\text{sas_erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
+        sas_erf\ $(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
         and complementary error function
         $\text{sas_erfc}(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
+        sas_erfc\ $(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
         The standard math functions erf(x) and erfc(x) are slower and broken
 …
     sas_J0(x):
         Bessel function of the first kind $\text{sas_J0}(x)=J_0(x)$ where
+        Bessel function of the first kind sas_J0\ $(x)=J_0(x)$ where
         $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$.
 …
     sas_J1(x):
         Bessel function of the first kind  $\text{sas_J1}(x)=J_1(x)$ where
+        Bessel function of the first kind  sas_J1\ $(x)=J_1(x)$ where
         $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$.
 …
     sas_JN(n, x):
         Bessel function of the first kind and integer order $n$:
         $\text{sas_JN}(n, x)=J_n(x)$ where
+        Bessel function of the first kind and integer order $n$,
+        sas_JN\ $(n, x) =J_n(x)$ where
         $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$.
         If $n$ = 0 or 1, it uses sas_J0(x) or sas_J1(x), respectively.
+        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively.
         The standard math function jn(n, x) is not available on all platforms.
 …
     sas_Si(x):
         Sine integral $\text{Si}(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
+        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
         This function uses Taylor series for small and large arguments:
 …
     sas_3j1x_x(x):
         Spherical Bessel form
         $\text{sph_j1c}(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
+        sph_j1c\ $(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
         with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical
         Bessel function of the first kind and first order.
 …
     sas_2J1x_x(x):
         Bessel form $\text{sas_J1c}(x) = 2 J_1(x)/x$, with a limiting value
+        Bessel form sas_J1c\ $(x) = 2 J_1(x)/x$, with a limiting value
         of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind
         and first order.
 …
 across multiple parameters can be very slow).
 If your model has 2D orientational calculation, then you should also
+If your model has 2D orientation calculation, then you should also
 test with::

doc/guide/resolution.rst

-                      rf8a2baa
+                      r1f058ea
 resolution contribution into a model calculation/simulation (which by definition
 will be exact) to make it more representative of what has been measured
 experimentally - a process called *smearing*. sasmodels does the latter.
+experimentally - a process called *smearing*. Sasmodels does the latter.
 Both smearing and desmearing rely on functions to describe the resolution
 effect. sasmodels provides three smearing algorithms:
+effect. Sasmodels provides three smearing algorithms:
 *  *Slit Smearing*
 …
 For discrete $q$ values, at the $q$ values of the data points and at the $q$
 values extended up to $q_N = q_i + \Delta q_v$ the smeared
+values extended up to $q_N = q_i + \Delta q_u$ the smeared
 intensity can be approximately calculated as
 …
 elliptical Gaussian distribution. The $A$ is a normalization factor.
 .. figure:: resolution_2d_rotation.gif
+.. figure:: resolution_2d_rotation.png
     Coordinate axis rotation for 2D resolution calculation.

doc/guide/scripting.rst

-                      r2e66ef5
+                      r4aa5dce
     $ bumps example/model.py --preview
+Note that bumps and sasmodels are included as part of the SasView
+distribution.  On windows, bumps can be called from the cmd prompt
+as follows::
+    SasViewCom bumps.cli example/model.py --preview
 Using sasmodels directly
 ========================
 …
     plt.loglog(q, Iq)
     plt.show()
+On windows, this can be called from the cmd prompt using sasview as::
+    SasViewCom example/cylinder_eval.py

doc/rst_prolog

-                      ra0fb06a
+                      r30b60d2
 .. Set up some substitutions to make life easier...
-.. Remove |biggamma|, etc. when they are no longer needed.
-.. |alpha| unicode:: U+03B1
-.. |beta| unicode:: U+03B2
-.. |gamma| unicode:: U+03B3
-.. |delta| unicode:: U+03B4
-.. |epsilon| unicode:: U+03B5
-.. |zeta| unicode:: U+03B6
-.. |eta| unicode:: U+03B7
-.. |theta| unicode:: U+03B8
-.. |iota| unicode:: U+03B9
-.. |kappa| unicode:: U+03BA
-.. |lambda| unicode:: U+03BB
-.. |mu| unicode:: U+03BC
-.. |nu| unicode:: U+03BD
-.. |xi| unicode:: U+03BE
-.. |omicron| unicode:: U+03BF
-.. |pi| unicode:: U+03C0
-.. |rho| unicode:: U+03C1
-.. |sigma| unicode:: U+03C3
-.. |tau| unicode:: U+03C4
-.. |upsilon| unicode:: U+03C5
-.. |phi| unicode:: U+03C6
-.. |chi| unicode:: U+03C7
-.. |psi| unicode:: U+03C8
-.. |omega| unicode:: U+03C9
-.. |biggamma| unicode:: U+0393
-.. |bigdelta| unicode:: U+0394
-.. |bigzeta| unicode:: U+039E
-.. |bigpsi| unicode:: U+03A8
-.. |bigphi| unicode:: U+03A6
-.. |bigsigma| unicode:: U+03A3
-.. |Gamma| unicode:: U+0393
-.. |Delta| unicode:: U+0394
-.. |Zeta| unicode:: U+039E
-.. |Psi| unicode:: U+03A8
-.. |drho| replace:: |Delta|\ |rho|
 .. |Ang| unicode:: U+212B
 .. |Ang^-1| replace:: |Ang|\ :sup:`-1`
 …
 .. |cm^-3| replace:: cm\ :sup:`-3`
 .. |sr^-1| replace:: sr\ :sup:`-1`
-.. |P0| replace:: P\ :sub:`0`\
-.. |A2| replace:: A\ :sub:`2`\
-.. |equiv| unicode:: U+2261
-.. |noteql| unicode:: U+2260
-.. |TM| unicode:: U+2122
 .. |cdot| unicode:: U+00B7

example/oriented_usans.py

r1cd24b4	r74b0495
6	6
7	7	# Spherical particle data, not ellipsoids
8		sans, usans = load_data('~~../../sasview/sasview/test/1d_data/latex_smeared.xm~~l')
	8	sans, usans = load_data('latex_smeared.xml', index='all')
9	9	usans.qmin, usans.qmax = np.min(usans.x), np.max(usans.x)
10	10	usans.mask = (usans.x < 0.0)

example/simul_fit.py

-                      r1a4d4c0
+                      r74b0495
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-# To Sasview/documents/scripts
 from bumps.names import *
 from sasmodels.core import load_model
 …
 from sasmodels.data import load_data, plot_data
+# latex data, same sample usans and sans
+# particles radius ~2300, uniform dispersity
+datasets = load_data('latex_smeared.xml', index='all')
+#[print(data) for data in datasets]
+""" IMPORT THE DATA USED """
+datafiles = ['latex_smeared_out_0.txt', 'latex_smeared_out_1.txt']
+datasets = [load_data(el) for el in datafiles]
+# A single sphere model to share between the datasets.  We will use
+# FreeVariables below to set the parameters that are independent between
+# the datasets.
+kernel = load_model('sphere')
+pars = dict(scale=0.01, background=0.0, sld=5.0, sld_solvent=0.0, radius=1500.,
+            #radius_pd=0.1, radius_pd_n=35,
+            )
+model = Model(kernel, **pars)
+for data in datasets:
+    data.qmin = 0.0
+    data.qmax = 10.0
+# radius and polydispersity (if any) are shared
+model.radius.range(0, inf)
+#model.radius_pd.range(0, 1)
+#sphere model
+kernel = load_model('sphere', dtype="single")
+pars = dict(scale=0.01, background=0.0, sld=1.0, sld_solvent=6.0, radius=1500.)
+model = Model(kernel, **pars)
+model.radius.range(0, inf)
+#model.background.range(-inf, inf)
+#model.scale.range(0, inf)
+model.sld.range(-inf, inf)
+model.sld_solvent.range(-inf, inf)
+# Contrast and dilution are the same for both measurements, but are not
+# separable with a single measurement (i.e., I(q) ~ F(q) contrast^2 Vf),
+# so fit one of scale, sld or solvent sld.  With absolute scaling from
+# data reduction, can use the same parameter for both datasets.
+model.scale.range(0, inf)
+#model.sld.range(-inf, inf)
+#model.sld_solvent.range(-inf, inf)
+# Background is different for sans and usans so set it as a free variable
+# in the model.
 free = FreeVariables(
     names=[data.filename for data in datasets],
+    names=[data.run[0] for data in datasets],
     background=model.background,
-    scale=model.scale,
+    )
 free.background.range(-inf, inf)
-free.scale.range(0, inf)
+M = [Experiment(data=data, model=model) for data in datasets]
+# Note: can access the parameters for the individual models using
+# free.background[0] and free.background[1], setting constraints or
+# ranges as appropriate.
+# For more complex systems where different datasets require independent models,
+# separate models can be defined, with parameters tied together using
+# constraint expressions.  For example, the following could be used to fit
+# data set 1 to spheres and data set 2 to cylinders of the same volume:
+#    model1 = Model(load_model('sphere'))
+#    model2 = Model(load_model('cylinder'))
+#    model1.sld = model2.sld
+#    model1.sld_solvent = model2.sld_solvent
+#    model1.scale = model2.scale
+#    # set cylinders and spheres to the same volume
+#    model1.radius = (3/4*model2.radius**2*model2.length)**(1/3)
+#    model1.background.range(0, 2)
+#    model2.background.range(0, 2)
+# Setup the experiments, sharing the same model across all datasets.
+M = [Experiment(data=data, model=model, name=data.run[0]) for data in datasets]
 problem = FitProblem(M, freevars=free)
-print(problem._parameters)

explore/precision.py

-                      r487e695
+                      r237c9cf
         elif xrange == "linear":
             lin_min, lin_max, lin_steps = 1, 1000, 2000
+            lin_min, lin_max, lin_steps = 0.001, 2, 2000
         elif xrange == "log":
             log_min, log_max, log_steps = -3, 5, 400
 …
     np_function=lambda x: np.fmod(x, 2*np.pi),
     ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
+)
+add_function(
+    name="debye",
+    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
+        const double x = qsq;
+        if (0) { // 0.36 single
+            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}]
+            return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1);
+        } else if (0) { // 1.0 for single
+            // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}]
+            const double A1=1./24., A2=1./84, A3=-1./3360;
+            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
+            // 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.;
+            const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.;
+            return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.)
+                  /((((B4*x + B3)*x + B2)*x + B1)*x + 1.);
+        } else { // 1.0 for single, 0.5 for double
+            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
+            const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.;
+            const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.;
+            return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.)
+                  /(((((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
+        const double x = qsq;
+        const double C0 = +1.;
+        const double C1 = -1./3.;
+        const double C2 = +1./12.;
+        const double C3 = -1./60.;
+        const double C4 = +1./360.;
+        const double C5 = -1./2520.;
+        const double C6 = +1./20160.;
+        const double C7 = -1./181440.;
+        //return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
+        //return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
+        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
+    } else {
+        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq);
+    }
+    """, "sas_debye"),
+)

sasmodels/bumps_model.py

-                      r3330bb4
+                      r74b0495
     """
     _cache = None # type: Dict[str, np.ndarray]
     def __init__(self, data, model, cutoff=1e-5):
+    def __init__(self, data, model, cutoff=1e-5, name=None):
         # type: (Data, Model, float) -> None
         # remember inputs so we can inspect from outside
+        self.name = data.filename if name is None else name
         self.model = model
         self.cutoff = cutoff
 …
         """
         data, theory, resid = self._data, self.theory(), self.residuals()
+        plot_theory(data, theory, resid, view, Iq_calc=self.Iq_calc)
+        # TODO: hack to display oriented usans 2-D pattern
+        Iq_calc = self.Iq_calc if isinstance(self.Iq_calc, tuple) else None
+        plot_theory(data, theory, resid, view, Iq_calc=Iq_calc)
     def simulate_data(self, noise=None):

sasmodels/compare.py

-                      r630156b
+                      rced5bd2
 USAGE = """
+usage: compare.py model N1 N2 [options...] [key=val]
+Compare the speed and value for a model between the SasView original and the
+sasmodels rewrite.
+usage: sascomp model [options...] [key=val]
+Generate and compare SAS models.  If a single model is specified it shows
+a plot of that model.  Different models can be compared, or the same model
+with different parameters.  The same model with the same parameters can
+be compared with different calculation engines to see the effects of precision
+on the resultant values.
 model or model1,model2 are the names of the models to compare (see below).
-N1 is the number of times to run sasmodels (default=1).
-N2 is the number times to run sasview (default=1).
 Options (* for default):
+    -plot*/-noplot plots or suppress the plot of the model
+    === data generation ===
+    -data="path" uses q, dq from the data file
+    -noise=0 sets the measurement error dI/I
+    -res=0 sets the resolution width dQ/Q if calculating with resolution
     -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
+    -q=min:max alternative specification of qrange
     -nq=128 sets the number of Q points in the data set
+    -1d*/-2d computes 1d or 2d data
     -zero indicates that q=0 should be included
+    -1d*/-2d computes 1d or 2d data
+    === model parameters ===
     -preset*/-random[=seed] preset or random parameters
+    -sets=n generates n random datasets with the seed given by -random=seed
+    -pars/-nopars* prints the parameter set or not
+    -default/-demo* use demo vs default parameters
+    === calculation options ===
     -mono*/-poly force monodisperse or allow polydisperse demo parameters
+    -cutoff=1e-5* cutoff value for including a point in polydispersity
     -magnetic/-nonmagnetic* suppress magnetism
-    -cutoff=1e-5* cutoff value for including a point in polydispersity
-    -pars/-nopars* prints the parameter set or not
-    -abs/-rel* plot relative or absolute error
-    -linear/-log*/-q4 intensity scaling
-    -hist/-nohist* plot histogram of relative error
-    -res=0 sets the resolution width dQ/Q if calculating with resolution
     -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
+    -edit starts the parameter explorer
+    -default/-demo* use demo vs default parameters
+    -help/-html shows the model docs instead of running the model
+    -title="note" adds note to the plot title, after the model name
+    -data="path" uses q, dq from the data file
+Any two calculation engines can be selected for comparison:
+    -neval=1 sets the number of evals for more accurate timing
+    === precision options ===
+    -calc=default uses the default calcution precision
     -single/-double/-half/-fast sets an OpenCL calculation engine
     -single!/-double!/-quad! sets an OpenMP calculation engine
     -sasview sets the sasview calculation engine
+The default is -single -double.  Note that the interpretation of quad
+precision depends on architecture, and may vary from 64-bit to 128-bit,
+with 80-bit floats being common (1e-19 precision).
+    === plotting ===
+    -plot*/-noplot plots or suppress the plot of the model
+    -linear/-log*/-q4 intensity scaling on plots
+    -hist/-nohist* plot histogram of relative error
+    -abs/-rel* plot relative or absolute error
+    -title="note" adds note to the plot title, after the model name
+    === output options ===
+    -edit starts the parameter explorer
+    -help/-html shows the model docs instead of running the model
+The interpretation of quad precision depends on architecture, and may
+vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision).
+On unix and mac you may need single quotes around the DLL computation
+engines, such as -calc='single!,double!' since !, is treated as a history
+expansion request in the shell.
 Key=value pairs allow you to set specific values for the model parameters.
+Key=value1,value2 to compare different values of the same parameter.
+value can be an expression including other parameters
+Key=value1,value2 to compare different values of the same parameter. The
+value can be an expression including other parameters.
+Items later on the command line override those that appear earlier.
+Examples:
+    # compare single and double precision calculation for a barbell
+    sascomp barbell -calc=single,double
+    # generate 10 random lorentz models, with seed=27
+    sascomp lorentz -sets=10 -seed=27
+    # compare ellipsoid with R = R_polar = R_equatorial to sphere of radius R
+    sascomp sphere,ellipsoid radius_polar=radius radius_equatorial=radius
+    # model timing test requires multiple evals to perform the estimate
+    sascomp pringle -calc=single,double -timing=100,100 -noplot
 """
 …
 -------------------
+"""
+           + USAGE)
+""" + USAGE)
 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
 …
 def _randomize_one(model_info, p, v):
+def _randomize_one(model_info, name, value):
     # type: (ModelInfo, str, float) -> float
     # type: (ModelInfo, str, str) -> str
 …
     Randomize a single parameter.
     """
+    if any(p.endswith(s) for s in ('_pd', '_pd_n', '_pd_nsigma', '_pd_type')):
+        return v
+    # Find the parameter definition
+    for par in model_info.parameters.call_parameters:
+        if par.name == p:
+            break
+    else:
+        raise ValueError("unknown parameter %r"%p)
+    # Set the amount of polydispersity/angular dispersion, but by default pd_n
+    # is zero so there is no polydispersity.  This allows us to turn on/off
+    # pd by setting pd_n, and still have randomly generated values
+    if name.endswith('_pd'):
+        par = model_info.parameters[name[:-3]]
+        if par.type == 'orientation':
+            # Let oriention variation peak around 13 degrees; 95% < 42 degrees
+            return 180*np.random.beta(2.5, 20)
+        else:
+            # Let polydispersity peak around 15%; 95% < 0.4; max=100%
+            return np.random.beta(1.5, 7)
+    # pd is selected globally rather than per parameter, so set to 0 for no pd
+    # In particular, when multiple pd dimensions, want to decrease the number
+    # of points per dimension for faster computation
+    if name.endswith('_pd_n'):
+        return 0
+    # Don't mess with distribution type for now
+    if name.endswith('_pd_type'):
+        return 'gaussian'
+    # type-dependent value of number of sigmas; for gaussian use 3.
+    if name.endswith('_pd_nsigma'):
+        return 3.
+    # background in the range [0.01, 1]
+    if name == 'background':
+        return 10**np.random.uniform(-2, 0)
+    # scale defaults to 0.1% to 30% volume fraction
+    if name == 'scale':
+        return 10**np.random.uniform(-3, -0.5)
+    # 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))
+    limits = par.limits
+    if not np.isfinite(limits).all():
+        low, high = parameter_range(p, v)
+        limits = (max(limits[0], low), min(limits[1], high))
+    # If it is a fixed range, pick from it with equal probability.
+    # For logarithmic ranges, the model will have to override.
+    if np.isfinite(par.limits).all():
+        return np.random.uniform(*par.limits)
+    # If the paramter is marked as an sld use the range of neutron slds
+    # TODO: ought to randomly contrast match a pair of SLDs
+    if par.type == 'sld':
+        return np.random.uniform(-0.5, 12)
+    # Limit magnetic SLDs to a smaller range, from zero to iron=5/A^2
+    if par.name.startswith('M0:'):
+        return np.random.uniform(0, 5)
+    # Guess at the random length/radius/thickness.  In practice, all models
+    # are going to set their own reasonable ranges.
+    if par.type == 'volume':
+        if ('length' in par.name or
+                'radius' in par.name or
+                'thick' in par.name):
+            return 10**np.random.uniform(2, 4)
+    # In the absence of any other info, select a value in [0, 2v], or
+    # [-2|v|, 2|v|] if v is negative, or [0, 1] if v is zero.  Mostly the
+    # model random parameter generators will override this default.
+    low, high = parameter_range(par.name, value)
+    limits = (max(par.limits[0], low), min(par.limits[1], high))
     return np.random.uniform(*limits)
+def randomize_pars(model_info, pars, seed=None):
+    # type: (ModelInfo, ParameterSet, int) -> ParameterSet
+def _random_pd(model_info, pars):
+    pd = [p for p in model_info.parameters.kernel_parameters if p.polydisperse]
+    pd_volume = []
+    pd_oriented = []
+    for p in pd:
+        if p.type == 'orientation':
+            pd_oriented.append(p.name)
+        elif p.length_control is not None:
+            n = int(pars.get(p.length_control, 1) + 0.5)
+            pd_volume.extend(p.name+str(k+1) for k in range(n))
+        elif p.length > 1:
+            pd_volume.extend(p.name+str(k+1) for k in range(p.length))
+        else:
+            pd_volume.append(p.name)
+    u = np.random.rand()
+    n = len(pd_volume)
+    if u < 0.01 or n < 1:
+        pass  # 1% chance of no polydispersity
+    elif u < 0.86 or n < 2:
+        pars[np.random.choice(pd_volume)+"_pd_n"] = 35
+    elif u < 0.99 or n < 3:
+        choices = np.random.choice(len(pd_volume), size=2)
+        pars[pd_volume[choices[0]]+"_pd_n"] = 25
+        pars[pd_volume[choices[1]]+"_pd_n"] = 10
+    else:
+        choices = np.random.choice(len(pd_volume), size=3)
+        pars[pd_volume[choices[0]]+"_pd_n"] = 25
+        pars[pd_volume[choices[1]]+"_pd_n"] = 10
+        pars[pd_volume[choices[2]]+"_pd_n"] = 5
+    if pd_oriented:
+        pars['theta_pd_n'] = 20
+        if np.random.rand() < 0.1:
+            pars['phi_pd_n'] = 5
+        if np.random.rand() < 0.1:
+            if any(p.name == 'psi' for p in model_info.parameters.kernel_parameters):
+                #print("generating psi_pd_n")
+                pars['psi_pd_n'] = 5
+    ## Show selected polydispersity
+    #for name, value in pars.items():
+    #    if name.endswith('_pd_n') and value > 0:
+    #        print(name, value, pars.get(name[:-5], 0), pars.get(name[:-2], 0))
+def randomize_pars(model_info, pars):
+    # type: (ModelInfo, ParameterSet) -> ParameterSet
     """
     Generate random values for all of the parameters.
 …
     :func:`constrain_pars` needs to be called afterward..
     """
+    with push_seed(seed):
+        # Note: the sort guarantees order `of calls to random number generator
+        random_pars = dict((p, _randomize_one(model_info, p, v))
+                           for p, v in sorted(pars.items()))
+    # Note: the sort guarantees order of calls to random number generator
+    random_pars = dict((p, _randomize_one(model_info, p, v))
+                       for p, v in sorted(pars.items()))
+    if model_info.random is not None:
+        random_pars.update(model_info.random())
+    _random_pd(model_info, random_pars)
     return random_pars
 def constrain_pars(model_info, pars):
 …
     Warning: this updates the *pars* dictionary in place.
     """
+    # TODO: move the model specific code to the individual models
     name = model_info.id
     # if it is a product model, then just look at the form factor since
 …
     elif name == 'guinier':
         # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
+        # I(q) = A e^-(Rg^2 q^2/3) > e^-(30 ln 10)
+        # => ln A - (Rg^2 q^2/3) > -30 ln 10
+        # => Rg^2 q^2/3 < 30 ln 10 + ln A
+        # => Rg < sqrt(90 ln 10 + 3 ln A)/q
         #q_max = 0.2  # mid q maximum
         q_max = 1.0  # high q maximum
 …
     parameters = model_info.parameters
     magnetic = False
+    magnetic_pars = []
     for p in parameters.user_parameters(pars, is2d):
         if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')):
             continue
+        if p.id.startswith('up:') and not magnetic:
+        if p.id.startswith('up:'):
+            magnetic_pars.append("%s=%s"%(p.id, pars.get(p.id, p.default)))
             continue
         fields = dict(
 …
         lines.append(_format_par(p.name, **fields))
         magnetic = magnetic or fields['M0'] != 0.
+    if magnetic and magnetic_pars:
+        lines.append(" ".join(magnetic_pars))
     return "\n".join(lines)
 …
                 % (pd, n, nsigma, nsigma, pdtype)
     if M0 != 0.:
         line += "  M0:%.3f  mphi:%.1f  mtheta:%.1f" % (M0, mphi, mtheta)
+        line += "  M0:%.3f  mtheta:%.1f  mphi:%.1f" % (M0, mtheta, mphi)
     return line
 def suppress_pd(pars):
+def suppress_pd(pars, suppress=True):
     # type: (ParameterSet) -> ParameterSet
     """
+    Suppress theta_pd for now until the normalization is resolved.
+    May also suppress complete polydispersity of the model to test
+    models more quickly.
+    If suppress is True complete eliminate polydispersity of the model to test
+    models more quickly.  If suppress is False, make sure at least one
+    parameter is polydisperse, setting the first polydispersity parameter to
+% if no polydispersity is given (with no explicit demo parameters given
+    in the model, there will be no default polydispersity).
     """
     pars = pars.copy()
+    for p in pars:
+        if p.endswith("_pd_n"): pars[p] = 0
+    #print("pars=", pars)
+    if suppress:
+        for p in pars:
+            if p.endswith("_pd_n"):
+                pars[p] = 0
+    else:
+        any_pd = False
+        first_pd = None
+        for p in pars:
+            if p.endswith("_pd_n"):
+                pd = pars.get(p[:-2], 0.)
+                any_pd |= (pars[p] != 0 and pd != 0.)
+                if first_pd is None:
+                    first_pd = p
+        if not any_pd and first_pd is not None:
+            if pars[first_pd] == 0:
+                pars[first_pd] = 35
+            if first_pd[:-2] not in pars or pars[first_pd[:-2]] == 0:
+                pars[first_pd[:-2]] = 0.15
     return pars
 def suppress_magnetism(pars):
+def suppress_magnetism(pars, suppress=True):
     # type: (ParameterSet) -> ParameterSet
     """
+    Suppress theta_pd for now until the normalization is resolved.
+    May also suppress complete polydispersity of the model to test
+    models more quickly.
+    If suppress is True complete eliminate magnetism of the model to test
+    models more quickly.  If suppress is False, make sure at least one sld
+    parameter is magnetic, setting the first parameter to have a strong
+    magnetic sld (8/A^2) at 60 degrees (with no explicit demo parameters given
+    in the model, there will be no default magnetism).
     """
     pars = pars.copy()
+    for p in pars:
+        if p.startswith("M0:"): pars[p] = 0
+    if suppress:
+        for p in pars:
+            if p.startswith("M0:"):
+                pars[p] = 0
+    else:
+        any_mag = False
+        first_mag = None
+        for p in pars:
+            if p.startswith("M0:"):
+                any_mag |= (pars[p] != 0)
+                if first_mag is None:
+                    first_mag = p
+        if not any_mag and first_mag is not None:
+            pars[first_mag] = 8.
     return pars
 …
     model = core.build_model(model_info, dtype=dtype, platform="ocl")
     calculator = DirectModel(data, model, cutoff=cutoff)
     calculator.engine = "OCL%s"%DTYPE_MAP[dtype]
+    calculator.engine = "OCL%s"%DTYPE_MAP[str(model.dtype)]
     return calculator
 …
     model = core.build_model(model_info, dtype=dtype, platform="dll")
     calculator = DirectModel(data, model, cutoff=cutoff)
     calculator.engine = "OMP%s"%DTYPE_MAP[dtype]
+    calculator.engine = "OMP%s"%DTYPE_MAP[str(model.dtype)]
     return calculator
 …
     'accuracy', 'is2d' and 'view' parsed from the command line.
     """
     qmax, nq, res = opts['qmax'], opts['nq'], opts['res']
+    qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res']
     if opts['is2d']:
         q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray
 …
     else:
         if opts['view'] == 'log' and not opts['zero']:
+            qmax = math.log10(qmax)
+            q = np.logspace(qmax-3, qmax, nq)
+            q = np.logspace(math.log10(qmin), math.log10(qmax), nq)
         else:
             q = np.linspace(0.001*qmax, qmax, nq)
+            q = np.linspace(qmin, qmax, nq)
         if opts['zero']:
             q = np.hstack((0, q))
 …
     if dtype == 'sasview':
         return eval_sasview(model_info, data)
+    elif dtype.endswith('!'):
+    elif dtype is None or not dtype.endswith('!'):
+        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
+    else:
         return eval_ctypes(model_info, data, dtype=dtype[:-1], cutoff=cutoff)
-    else:
-        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
 def _show_invalid(data, theory):
 …
     parameters.
     """
+    result = run_models(opts, verbose=True)
+    if opts['plot']:  # Note: never called from explore
+        plot_models(opts, result, limits=limits)
+    limits = np.Inf, -np.Inf
+    for k in range(opts['sets']):
+        opts['pars'] = parse_pars(opts)
+        if opts['pars'] is None:
+            return
+        result = run_models(opts, verbose=True)
+        if opts['plot']:
+            limits = plot_models(opts, result, limits=limits, setnum=k)
+    if opts['plot']:
+        import matplotlib.pyplot as plt
+        plt.show()
 def run_models(opts, verbose=False):
     # type: (Dict[str, Any]) -> Dict[str, Any]
+    n_base, n_comp = opts['count']
+    pars, pars2 = opts['pars']
+    base, comp = opts['engines']
+    base_n, comp_n = opts['count']
+    base_pars, comp_pars = opts['pars']
     data = opts['data']
+    # silence the linter
+    base = opts['engines'][0] if n_base else None
+    comp = opts['engines'][1] if n_comp else None
+    comparison = comp is not None
     base_time = comp_time = None
 …
     # Base calculation
+    if n_base > 0:
+    try:
+        base_raw, base_time = time_calculation(base, base_pars, base_n)
+        base_value = np.ma.masked_invalid(base_raw)
+        if verbose:
+            print("%s t=%.2f ms, intensity=%.0f"
+                  % (base.engine, base_time, base_value.sum()))
+        _show_invalid(data, base_value)
+    except ImportError:
+        traceback.print_exc()
+    # Comparison calculation
+    if comparison:
         try:
+            base_raw, base_time = time_calculation(base, pars, n_base)
+            base_value = np.ma.masked_invalid(base_raw)
+            if verbose:
+                print("%s t=%.2f ms, intensity=%.0f"
+                      % (base.engine, base_time, base_value.sum()))
+            _show_invalid(data, base_value)
+        except ImportError:
+            traceback.print_exc()
+            n_base = 0
+    # Comparison calculation
+    if n_comp > 0:
+        try:
+            comp_raw, comp_time = time_calculation(comp, pars2, n_comp)
+            comp_raw, comp_time = time_calculation(comp, comp_pars, comp_n)
             comp_value = np.ma.masked_invalid(comp_raw)
             if verbose:
 …
         except ImportError:
             traceback.print_exc()
-            n_comp = 0
     # Compare, but only if computing both forms
     if n_base > 0 and n_comp > 0:
+    if comparison:
         resid = (base_value - comp_value)
         relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0)
 …
 def plot_models(opts, result, limits=None):
+def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0):
     # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
     base_value, comp_value= result['base_value'], result['comp_value']
+    base_value, comp_value = result['base_value'], result['comp_value']
     base_time, comp_time = result['base_time'], result['comp_time']
     resid, relerr = result['resid'], result['relerr']
     have_base, have_comp = (base_value is not None), (comp_value is not None)
+    base = opts['engines'][0] if have_base else None
+    comp = opts['engines'][1] if have_comp else None
+    base, comp = opts['engines']
     data = opts['data']
     use_data = (opts['datafile'] is not None) and (have_base ^ have_comp)
 …
     view = opts['view']
     import matplotlib.pyplot as plt
+    if limits is None and not use_data:
+        vmin, vmax = np.Inf, -np.Inf
+        if have_base:
+            vmin = min(vmin, base_value.min())
+            vmax = max(vmax, base_value.max())
+    vmin, vmax = limits
+    if have_base:
+        vmin = min(vmin, base_value.min())
+        vmax = max(vmax, base_value.max())
+    if have_comp:
+        vmin = min(vmin, comp_value.min())
+        vmax = max(vmax, comp_value.max())
+    limits = vmin, vmax
+    if have_base:
         if have_comp:
+            vmin = min(vmin, comp_value.min())
+            vmax = max(vmax, comp_value.max())
+        limits = vmin, vmax
+    if have_base:
+        if have_comp: plt.subplot(131)
+            plt.subplot(131)
         plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
         plt.title("%s t=%.2f ms"%(base.engine, base_time))
         #cbar_title = "log I"
     if have_comp:
+        if have_base: plt.subplot(132)
+        if have_base:
+            plt.subplot(132)
         if not opts['is2d'] and have_base:
             plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
 …
         else:
             err, errstr, errview = abs(relerr), "rel err", "log"
+            if (err == 0.).all():
+                errview = 'linear'
         if 0:  # 95% cutoff
             sorted = np.sort(err.flatten())
             cutoff = sorted[int(sorted.size*0.95)]
             err[err>cutoff] = cutoff
+            err[err > cutoff] = cutoff
         #err,errstr = base/comp,"ratio"
+        plot_theory(data, None, resid=err, view=errview, use_data=use_data)
+        if view == 'linear':
+            plt.xscale('linear')
+        plot_theory(data, None, resid=err, view=view, use_data=use_data)
+        plt.yscale(errview)
         plt.title("max %s = %.3g"%(errstr, abs(err).max()))
         #cbar_title = errstr if errview=="linear" else "log "+errstr
 …
         plt.title('Distribution of relative error between calculation engines')
-    if not opts['explore']:
-        plt.show()
     return limits
 # ===========================================================================
+#
+NAME_OPTIONS = set([
+# Set of command line options.
+# Normal options such as -plot/-noplot are specified as 'name'.
+# For options such as -nq=500 which require a value use 'name='.
+#
+OPTIONS = [
+    # Plotting
     'plot', 'noplot',
+    'half', 'fast', 'single', 'double',
+    'single!', 'double!', 'quad!', 'sasview',
+    'linear', 'log', 'q4',
+    'rel', 'abs',
+    'hist', 'nohist',
+    'title=',
+    # Data generation
+    'data=', 'noise=', 'res=', 'nq=', 'q=',
     'lowq', 'midq', 'highq', 'exq', 'zero',
     '2d', '1d',
+    'preset', 'random',
+    'poly', 'mono',
+    # Parameter set
+    'preset', 'random', 'random=', 'sets=',
+    'demo', 'default',  # TODO: remove demo/default
+    'nopars', 'pars',
+    # Calculation options
+    'poly', 'mono', 'cutoff=',
     'magnetic', 'nonmagnetic',
     'nopars', 'pars',
     'rel', 'abs',
+    'linear', 'log', 'q4',
     'hist', 'nohist',
     'edit', 'html', 'help',
     'demo', 'default',
     ])
+VALUE_OPTIONS = [
     # Note: random is both a name option and a value option
     'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data',
+    'accuracy=',
+    'neval=',  # for timing...
+    # Precision options
+    'calc=',
+    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
+    'sasview',  # TODO: remove sasview 3.x support
+    # Output options
+    'help', 'html', 'edit',
+    ]
+NAME_OPTIONS = set(k for k in OPTIONS if not k.endswith('='))
+VALUE_OPTIONS = [k[:-1] for k in OPTIONS if k.endswith('=')]
 def columnize(items, indent="", width=79):
 …
     # Plug in values given in demo
     if use_demo:
+    if use_demo and model_info.demo:
         pars.update(model_info.demo)
     return pars
 …
 # not sure about brackets [] {}
 # maybe one of the following @ ? ^ ! ,
 MODEL_SPLIT = ','
+PAR_SPLIT = ','
 def parse_opts(argv):
     # type: (List[str]) -> Dict[str, Any]
 …
               if not arg.startswith('-') and '=' in arg]
     positional_args = [arg for arg in argv
             if not arg.startswith('-') and '=' not in arg]
+                       if not arg.startswith('-') and '=' not in arg]
     models = "\n    ".join("%-15s"%v for v in MODELS)
     if len(positional_args) == 0:
 …
         print(columnize(MODELS, indent="  "))
         return None
-    if len(positional_args) > 3:
-        print("expected parameters: model N1 N2")
     invalid = [o[1:] for o in flags
 …
         return None
+    name = positional_args[0]
+    n1 = int(positional_args[1]) if len(positional_args) > 1 else 1
+    n2 = int(positional_args[2]) if len(positional_args) > 2 else 1
+    name = positional_args[-1]
     # pylint: disable=bad-whitespace
 …
         'view'      : 'log',
         'is2d'      : False,
+        'qmin'      : None,
         'qmax'      : 0.05,
         'nq'        : 128,
         'res'       : 0.0,
+        'noise'     : 0.0,
         'accuracy'  : 'Low',
         'cutoff'    : 0.0,
+        'cutoff'    : '0.0',
         'seed'      : -1,  # default to preset
         'mono'      : True,
 …
         'title'     : None,
         'datafile'  : None,
+        'sets'      : 0,
+        'engine'    : 'default',
+        'evals'     : '1',
+    }
-    engines = []
     for arg in flags:
         if arg == '-noplot':    opts['plot'] = False
 …
         elif arg == '-zero':    opts['zero'] = True
         elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
+        elif arg.startswith('-q='):
+            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')]
         elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
+        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
+        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
         elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
         elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:])
+        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:]
         elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
         elif arg.startswith('-title='):    opts['title'] = arg[7:]
         elif arg.startswith('-data='):     opts['datafile'] = arg[6:]
+        elif arg.startswith('-calc='):     opts['engine'] = arg[6:]
+        elif arg.startswith('-neval='):    opts['evals'] = arg[7:]
         elif arg == '-random':  opts['seed'] = np.random.randint(1000000)
         elif arg == '-preset':  opts['seed'] = -1
 …
         elif arg == '-rel':     opts['rel_err'] = True
         elif arg == '-abs':     opts['rel_err'] = False
         elif arg == '-half':    engines.append(arg[1:])
         elif arg == '-fast':    engines.append(arg[1:])
         elif arg == '-single':  engines.append(arg[1:])
         elif arg == '-double':  engines.append(arg[1:])
         elif arg == '-single!': engines.append(arg[1:])
         elif arg == '-double!': engines.append(arg[1:])
         elif arg == '-quad!':   engines.append(arg[1:])
         elif arg == '-sasview': engines.append(arg[1:])
+        elif arg == '-half':    opts['engine'] = 'half'
+        elif arg == '-fast':    opts['engine'] = 'fast'
+        elif arg == '-single':  opts['engine'] = 'single'
+        elif arg == '-double':  opts['engine'] = 'double'
+        elif arg == '-single!': opts['engine'] = 'single!'
+        elif arg == '-double!': opts['engine'] = 'double!'
+        elif arg == '-quad!':   opts['engine'] = 'quad!'
+        elif arg == '-sasview': opts['engine'] = 'sasview'
         elif arg == '-edit':    opts['explore'] = True
         elif arg == '-demo':    opts['use_demo'] = True
         elif arg == '-default':    opts['use_demo'] = False
+        elif arg == '-default': opts['use_demo'] = False
         elif arg == '-html':    opts['html'] = True
         elif arg == '-help':    opts['html'] = True
     # pylint: enable=bad-whitespace
+    if MODEL_SPLIT in name:
+        name, name2 = name.split(MODEL_SPLIT, 2)
+    # Magnetism forces 2D for now
+    if opts['magnetic']:
+        opts['is2d'] = True
+    # Force random if sets is used
+    if opts['sets'] >= 1 and opts['seed'] < 0:
+        opts['seed'] = np.random.randint(1000000)
+    if opts['sets'] == 0:
+        opts['sets'] = 1
+    # Create the computational engines
+    if opts['qmin'] is None:
+        opts['qmin'] = 0.001*opts['qmax']
+    if opts['datafile'] is not None:
+        data = load_data(os.path.expanduser(opts['datafile']))
     else:
+        name2 = name
+        data, _ = make_data(opts)
+    comparison = any(PAR_SPLIT in v for v in values)
+    if PAR_SPLIT in name:
+        names = name.split(PAR_SPLIT, 2)
+        comparison = True
+    else:
+        names = [name]*2
     try:
+        model_info = core.load_model_info(name)
+        model_info2 = core.load_model_info(name2) if name2 != name else model_info
+        model_info = [core.load_model_info(k) for k in names]
     except ImportError as exc:
         print(str(exc))
         print("Could not find model; use one of:\n    " + models)
         return None
+    if PAR_SPLIT in opts['engine']:
+        engine_types = opts['engine'].split(PAR_SPLIT, 2)
+        comparison = True
+    else:
+        engine_types = [opts['engine']]*2
+    if PAR_SPLIT in opts['evals']:
+        evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        evals = [int(opts['evals'])]*2
+    if PAR_SPLIT in opts['cutoff']:
+        cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        cutoff = [float(opts['cutoff'])]*2
+    base = make_engine(model_info[0], data, engine_types[0], cutoff[0])
+    if comparison:
+        comp = make_engine(model_info[1], data, engine_types[1], cutoff[1])
+    else:
+        comp = None
+    # pylint: disable=bad-whitespace
+    # Remember it all
+    opts.update({
+        'data'      : data,
+        'name'      : names,
+        'def'       : model_info,
+        'count'     : evals,
+        'engines'   : [base, comp],
+        'values'    : values,
+    })
+    # pylint: enable=bad-whitespace
+    return opts
+def parse_pars(opts):
+    model_info, model_info2 = opts['def']
     # Get demo parameters from model definition, or use default parameters
 …
     #pars.update(set_pars)  # set value before random to control range
     if opts['seed'] > -1:
         pars = randomize_pars(model_info, pars, seed=opts['seed'])
+        pars = randomize_pars(model_info, pars)
         if model_info != model_info2:
             pars2 = randomize_pars(model_info2, pars2, seed=opts['seed'])
+            pars2 = randomize_pars(model_info2, pars2)
             # Share values for parameters with the same name
             for k, v in pars.items():
 …
         constrain_pars(model_info, pars)
         constrain_pars(model_info2, pars2)
+        print("Randomize using -random=%i"%opts['seed'])
+    if opts['mono']:
+        pars = suppress_pd(pars)
+        pars2 = suppress_pd(pars2)
+    if not opts['magnetic']:
+        pars = suppress_magnetism(pars)
+        pars2 = suppress_magnetism(pars2)
+    pars = suppress_pd(pars, opts['mono'])
+    pars2 = suppress_pd(pars2, opts['mono'])
+    pars = suppress_magnetism(pars, not opts['magnetic'])
+    pars2 = suppress_magnetism(pars2, not opts['magnetic'])
     # Fill in parameters given on the command line
     presets = {}
     presets2 = {}
     for arg in values:
+    for arg in opts['values']:
         k, v = arg.split('=', 1)
         if k not in pars and k not in pars2:
 …
             print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
             return None
         v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v)
+        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v,v)
         if v1 and k in pars:
             presets[k] = float(v1) if isnumber(v1) else v1
 …
     context['np'] = np
     context.update(pars)
     context.update((k,v) for k,v in presets.items() if isinstance(v, float))
+    context.update((k, v) for k, v in presets.items() if isinstance(v, float))
     for k, v in presets.items():
         if not isinstance(v, float) and not k.endswith('_type'):
             presets[k] = eval(v, context)
     context.update(presets)
     context.update((k,v) for k,v in presets2.items() if isinstance(v, float))
+    context.update((k, v) for k, v in presets2.items() if isinstance(v, float))
     for k, v in presets2.items():
         if not isinstance(v, float) and not k.endswith('_type'):
 …
     #import pprint; pprint.pprint(model_info)
-    same_model = name == name2 and pars == pars
-    if len(engines) == 0:
-        if same_model:
-            engines.extend(['single', 'double'])
-        else:
-            engines.extend(['single', 'single'])
-    elif len(engines) == 1:
-        if not same_model:
-            engines.append(engines[0])
-        elif engines[0] == 'double':
-            engines.append('single')
-        else:
-            engines.append('double')
-    elif len(engines) > 2:
-        del engines[2:]
-    use_sasview = any(engine == 'sasview' and count > 0
-                      for engine, count in zip(engines, [n1, n2]))
-    if use_sasview:
-        constrain_new_to_old(model_info, pars)
-        constrain_new_to_old(model_info2, pars2)
     if opts['show_pars']:
         if not same_model:
+        if model_info.name != model_info2.name or pars != pars2:
             print("==== %s ====="%model_info.name)
             print(str(parlist(model_info, pars, opts['is2d'])))
 …
             print(str(parlist(model_info, pars, opts['is2d'])))
+    # Create the computational engines
+    if opts['datafile'] is not None:
+        data = load_data(os.path.expanduser(opts['datafile']))
+    else:
+        data, _ = make_data(opts)
+    if n1:
+        base = make_engine(model_info, data, engines[0], opts['cutoff'])
+    else:
+        base = None
+    if n2:
+        comp = make_engine(model_info2, data, engines[1], opts['cutoff'])
+    else:
+        comp = None
+    # pylint: disable=bad-whitespace
+    # Remember it all
+    opts.update({
+        'data'      : data,
+        'name'      : [name, name2],
+        'def'       : [model_info, model_info2],
+        'count'     : [n1, n2],
+        'presets'   : [presets, presets2],
+        'pars'      : [pars, pars2],
+        'engines'   : [base, comp],
+    })
+    # pylint: enable=bad-whitespace
+    return opts
+    return pars, pars2
 def show_docs(opts):
 …
     model = Explore(opts)
     problem = FitProblem(model)
+    frame = AppFrame(parent=None, title="explore", size=(1000,700))
+    if not is_mac: frame.Show()
+    frame = AppFrame(parent=None, title="explore", size=(1000, 700))
+    if not is_mac:
+        frame.Show()
     frame.panel.set_model(model=problem)
     frame.panel.Layout()
 …
     if is_mac: frame.Show()
     # If running withing an app, start the main loop
+    if app: app.MainLoop()
+    if app:
+        app.MainLoop()
 class Explore(object):
 …
         config_matplotlib()
         self.opts = opts
+        opts['pars'] = list(opts['pars'])
         p1, p2 = opts['pars']
         m1, m2 = opts['def']
 …
         self.starting_values = dict((k, v.value) for k, v in pars.items())
         self.pd_types = pd_types
         self.limits = None
+        self.limits = np.Inf, -np.Inf
     def revert_values(self):
 …
     opts = parse_opts(argv)
     if opts is not None:
+        if opts['seed'] > -1:
+            print("Randomize using -random=%i"%opts['seed'])
+            np.random.seed(opts['seed'])
         if opts['html']:
             show_docs(opts)
         elif opts['explore']:
+            opts['pars'] = parse_pars(opts)
+            if opts['pars'] is None:
+                return
             explore(opts)
         else:

sasmodels/conversion_table.py

-                      r0c2da4b
+                      r505d0ad
 account for that.
+Usage:
+<old_Sasview_version> : {
+    <new_model_name> : [
+        <old_model_name> ,
+        {
+            <new_param_name_1> : <old_param_name_1>,
+            ...
+            <new_param_name_n> : <old_param_name_n>
+        }
+    ]
+}
+Usage::
+    <old_Sasview_version> : {
+        <new_model_name> : [
+            <old_model_name> ,
+            {
+                <new_param_name_1> : <old_param_name_1>,
+                ...
+                <new_param_name_n> : <old_param_name_n>
+            }
+        ]
+    }
 Any future parameter and model name changes can and should be given in this

sasmodels/core.py

-                      r2e66ef5
+                      ra85a569
 import os
+import re
 from os.path import basename, dirname, join as joinpath
 from glob import glob
 …
 from . import kernelpy
 from . import kerneldll
+from . import custom
 if os.environ.get("SAS_OPENCL", "").lower() == "none":
 …
     except Exception:
         HAVE_OPENCL = False
+CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "")
+if CUSTOM_MODEL_PATH == "":
+    path = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models")
+    if not os.path.isdir(path):
+        os.makedirs(path)
+    CUSTOM_MODEL_PATH = path
 try:
 …
                        dtype=dtype, platform=platform)
+def load_model_info(model_name):
+def load_model_info(model_string):
     # type: (str) -> modelinfo.ModelInfo
     """
     Load a model definition given the model name.
+    *model_name* is the name of the model, or perhaps a model expression
+    such as sphere*hardsphere or sphere+cylinder.
+    *model_string* is the name of the model, or perhaps a model expression
+    such as sphere*cylinder or sphere+cylinder. Use '@' for a structure
+    factor product, e.g. sphere@hardsphere. Custom models can be specified by
+    prefixing the model name with 'custom.', e.g. 'custom.MyModel+sphere'.
     This returns a handle to the module defining the model.  This can be
     used with functions in generate to build the docs or extract model info.
     """
+    parts = model_name.split('+')
+    if len(parts) > 1:
+        model_info_list = [load_model_info(p) for p in parts]
+        return mixture.make_mixture_info(model_info_list)
+    parts = model_name.split('*')
+    if len(parts) > 1:
+        if len(parts) > 2:
+            raise ValueError("use P*S to apply structure factor S to model P")
+        P_info, Q_info = [load_model_info(p) for p in parts]
+    if '@' in model_string:
+        parts = model_string.split('@')
+        if len(parts) != 2:
+            raise ValueError("Use P@S to apply a structure factor S to model P")
+        P_info, Q_info = [load_model_info(part) for part in parts]
         return product.make_product_info(P_info, Q_info)
+    kernel_module = generate.load_kernel_module(model_name)
+    return modelinfo.make_model_info(kernel_module)
+    product_parts = []
+    addition_parts = []
+    addition_parts_names = model_string.split('+')
+    if len(addition_parts_names) >= 2:
+        addition_parts = [load_model_info(part) for part in addition_parts_names]
+    elif len(addition_parts_names) == 1:
+        product_parts_names = model_string.split('*')
+        if len(product_parts_names) >= 2:
+            product_parts = [load_model_info(part) for part in product_parts_names]
+        elif len(product_parts_names) == 1:
+            if "custom." in product_parts_names[0]:
+                # Extract ModelName from "custom.ModelName"
+                pattern = "custom.([A-Za-z0-9_-]+)"
+                result = re.match(pattern, product_parts_names[0])
+                if result is None:
+                    raise ValueError("Model name in invalid format: " + product_parts_names[0])
+                model_name = result.group(1)
+                # Use ModelName to find the path to the custom model file
+                model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py")
+                if not os.path.isfile(model_path):
+                    raise ValueError("The model file {} doesn't exist".format(model_path))
+                kernel_module = custom.load_custom_kernel_module(model_path)
+                return modelinfo.make_model_info(kernel_module)
+            # Model is a core model
+            kernel_module = generate.load_kernel_module(product_parts_names[0])
+            return modelinfo.make_model_info(kernel_module)
+    model = None
+    if len(product_parts) > 1:
+        model = mixture.make_mixture_info(product_parts, operation='*')
+    if len(addition_parts) > 1:
+        if model is not None:
+            addition_parts.append(model)
+        model = mixture.make_mixture_info(addition_parts, operation='+')
+    return model
 …
     # Convert dtype string to numpy dtype.
     if dtype is None:
+    if dtype is None or dtype == "default":
         numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single
                        else generate.F64)

sasmodels/custom/init.py

-                      r3c852d4
+                      r0f48f1e
         if fullname in sys.modules:
             del sys.modules[fullname]
+        if path.endswith(".py") and os.path.exists(path) and os.path.exists(path+"c"):
+            # remove automatic pyc file before loading a py file
+            os.unlink(path+"c")
         module = imp.load_source(fullname, os.path.expanduser(path))
-        #os.unlink(path+"c")  # remove the automatic pyc file
         return module

sasmodels/data.py

-                      r630156b
+                      rba7302a
     Data = Union["Data1D", "Data2D", "SesansData"]
 def load_data(filename):
+def load_data(filename, index=0):
     # type: (str) -> Data
     """
 …
         filename, indexstr = filename[:-1].split('[')
         index = int(indexstr)
-    else:
-        index = None
     datasets = loader.load(filename)
     if datasets is None:
+    if not datasets:  # None or []
         raise IOError("Data %r could not be loaded" % filename)
     if not isinstance(datasets, list):
         datasets = [datasets]
+    if index is None and len(datasets) > 1:
+        raise ValueError("Need to specify filename[index] for multipart data")
+    data = datasets[index if index is not None else 0]
+    if hasattr(data, 'x'):
+        data.qmin, data.qmax = data.x.min(), data.x.max()
+        data.mask = (np.isnan(data.y) if data.y is not None
+                     else np.zeros_like(data.x, dtype='bool'))
+    elif hasattr(data, 'qx_data'):
+        data.mask = ~data.mask
+    return data
+    for data in datasets:
+        if hasattr(data, 'x'):
+            data.qmin, data.qmax = data.x.min(), data.x.max()
+            data.mask = (np.isnan(data.y) if data.y is not None
+                        else np.zeros_like(data.x, dtype='bool'))
+        elif hasattr(data, 'qx_data'):
+            data.mask = ~data.mask
+    return datasets[index] if index != 'all' else datasets
 …
     scale = data.x**4 if view == 'q4' else 1.0
+    xscale = yscale = 'linear' if view == 'linear' else 'log'
     if use_data or use_theory:
 …
             plt.ylim(*limits)
+        plt.xscale('linear' if not some_present or non_positive_x
+                   else view if view is not None
+                   else 'log')
+        plt.yscale('linear'
+                   if view == 'q4' or not some_present or not all_positive
+                   else view if view is not None
+                   else 'log')
+        xscale = ('linear' if not some_present or non_positive_x
+                  else view if view is not None
+                  else 'log')
+        yscale = ('linear'
+                  if view == 'q4' or not some_present or not all_positive
+                  else view if view is not None
+                  else 'log')
+        plt.xscale(xscale)
         plt.xlabel("$q$/A$^{-1}$")
+        plt.yscale(yscale)
         plt.ylabel('$I(q)$')
         title = ("data and model" if use_theory and use_data
 …
         plt.xlabel("$q$/A$^{-1}$")
         plt.ylabel('residuals')
-        plt.xscale('linear')
         plt.title('(model - Iq)/dIq')
+        plt.xscale(xscale)
+        plt.yscale('linear')

sasmodels/direct_model.py

-                      ra769b54
+                      rd1ff3a5
         self.Iq = y
         if self.data_type in ('Iq', 'Iq-oriented'):
+            if self._data.y is None:
+                self._data.y = np.empty(len(self._data.x), 'd')
+            if self._data.dy is None:
+                self._data.dy = np.empty(len(self._data.x), 'd')
             self._data.dy[self.index] = dy
             self._data.y[self.index] = y
         elif self.data_type == 'Iqxy':
+            if self._data.data is None:
+                self._data.data = np.empty_like(self._data.qx_data, 'd')
+            if self._data.err_data is None:
+                self._data.err_data = np.empty_like(self._data.qx_data, 'd')
             self._data.data[self.index] = y
+            self._data.err_data[self.index] = dy
         elif self.data_type == 'sesans':
+            if self._data.y is None:
+                self._data.y = np.empty(len(self._data.x), 'd')
             self._data.y[self.index] = y
         else:
 …
         # TODO: extend plotting of calculate Iq to other measurement types
         # TODO: refactor so we don't store the result in the model
         self.Iq_calc = None
+        self.Iq_calc = Iq_calc
         if self.data_type == 'sesans':
             Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True)

sasmodels/generate.py

-                      rbb4b509
+                      r6db17bd
 # Conversion from units defined in the parameter table for each model
 # to units displayed in the sphinx documentation.
+# This section associates the unit with the macro to use to produce the LaTex
+# code.  The macro itself needs to be defined in sasmodels/doc/rst_prolog.
+#
+# NOTE: there is an RST_PROLOG at the end of this file which is NOT
+# used for the bundled documentation. Still as long as we are defining the macros
+# in two places any new addition should define the macro in both places.
 RST_UNITS = {
     "Ang": "|Ang|",
 …
     "1/Ang^2": "|Ang^-2|",
     "Ang^3": "|Ang^3|",
+    "Ang^2": "|Ang^2|",
     "1e15/cm^3": "|1e15cm^3|",
     "Ang^3/mol": "|Ang^3|/mol",
 …
     """
     # Note: need 0xffffffff&val to force an unsigned 32-bit number
+    try:
+        source = source.encode('utf8')
+    except AttributeError: # bytes has no encode attribute in python 3
+        pass
     return "%08X"%(0xffffffff&crc32(source))
 …
 .. |cm^-3| replace:: cm\ :sup:`-3`
 .. |sr^-1| replace:: sr\ :sup:`-1`
-.. |P0| replace:: P\ :sub:`0`\
-.. |equiv| unicode:: U+2261
-.. |noteql| unicode:: U+2260
-.. |TM| unicode:: U+2122
 .. |cdot| unicode:: U+00B7

sasmodels/kernelpy.py

-                      rdbfd471
+                      r3b32f8d
 from __future__ import division, print_function
+import logging
 import numpy as np  # type: ignore
 from numpy import pi, sin, cos  #type: ignore
 …
 try:
     from typing import Union, Callable
 except:
+except ImportError:
     pass
 else:
 …
         _create_default_functions(model_info)
         self.info = model_info
+        self.dtype = np.dtype('d')
+        logging.info("load python model " + self.info.name)
     def make_kernel(self, q_vectors):
 …
             # INVALID expression like the C models, but that is too expensive.
             Iq = np.asarray(form(), 'd')
+            if np.isnan(Iq).any(): continue
+            if np.isnan(Iq).any():
+                continue
             # update value and norm
 …
         default_Iqxy.vectorized = True
         model_info.Iqxy = default_Iqxy

sasmodels/list_pars.py

-                      r3330bb4
+                      r72be531
 import sys
+import argparse
 from .core import load_model_info, list_models
 from .compare import columnize
 def find_pars():
+def find_pars(type=None):
     """
     Find all parameters in all models.
 …
         model_info = load_model_info(name)
         for p in model_info.parameters.kernel_parameters:
+            partable.setdefault(p.name, [])
+            partable[p.name].append(name)
+            if type is None or p.type == type:
+                partable.setdefault(p.name, [])
+                partable[p.name].append(name)
     return partable
 def list_pars(names_only=True):
+def list_pars(names_only=True, type=None):
     """
     Print all parameters in all models.
 …
     occurs in.
     """
     partable = find_pars()
+    partable = find_pars(type)
     if names_only:
         print(columnize(list(sorted(partable.keys()))))
 …
     Program to list the parameters used across all models.
     """
+    if len(sys.argv) == 2 and sys.argv[1] == '-v':
+        verbose = True
+    elif len(sys.argv) == 1:
+        verbose = False
+    else:
+        print(__doc__)
+        sys.exit(1)
+    list_pars(names_only=not verbose)
+    parser = argparse.ArgumentParser(
+        description="Find all parameters in all models",
+        )
+    parser.add_argument(
+        '-v', '--verbose',
+        action='store_true',
+        help="list models which use this argument")
+    parser.add_argument(
+        'type', default="any", nargs='?',
+        metavar="volume|orientation|sld|none|any",
+        choices=['volume', 'orientation', 'sld', None, 'any'],
+        type=lambda v: None if v == 'any' else '' if v == 'none' else v,
+        help="only list arguments of the given type")
+    args = parser.parse_args()
+    list_pars(names_only=not args.verbose, type=args.type)
 if __name__ == "__main__":

sasmodels/mixture.py

-                      r6dc78e4
+                      r0edb6c1
     pass
 def make_mixture_info(parts):
+def make_mixture_info(parts, operation='+'):
     # type: (List[ModelInfo]) -> ModelInfo
     """
     Create info block for mixture model.
     """
-    flatten = []
-    for part in parts:
-        if part.composition and part.composition[0] == 'mixture':
-            flatten.extend(part.composition[1])
-        else:
-            flatten.append(part)
-    parts = flatten
     # Build new parameter list
     combined_pars = []
+    demo = {}
+    for k, part in enumerate(parts):
+    model_num = 0
+    all_parts = copy(parts)
+    is_flat = False
+    while not is_flat:
+        is_flat = True
+        for part in all_parts:
+            if part.composition and part.composition[0] == 'mixture' and \
+                len(part.composition[1]) > 1:
+                all_parts += part.composition[1]
+                all_parts.remove(part)
+                is_flat = False
+    # When creating a mixture model that is a sum of product models (ie (1*2)+(3*4))
+    # the parameters for models 1 & 2 will be prefixed with A & B respectively,
+    # but so will the parameters for models 3 & 4. We need to rename models 3 & 4
+    # so that they are prefixed with C & D to avoid overlap of parameter names.
+    used_prefixes = []
+    for part in parts:
+        i = 0
+        if part.composition and part.composition[0] == 'mixture':
+            npars_list = [info.parameters.npars for info in part.composition[1]]
+            for npars in npars_list:
+                # List of params of one of the constituent models of part
+                submodel_pars = part.parameters.kernel_parameters[i:i+npars]
+                # Prefix of the constituent model
+                prefix = submodel_pars[0].name[0]
+                if prefix not in used_prefixes: # Haven't seen this prefix so far
+                    used_prefixes.append(prefix)
+                    i += npars
+                    continue
+                while prefix in used_prefixes:
+                    # This prefix has been already used, so change it to the
+                    # next letter that hasn't been used
+                    prefix = chr(ord(prefix) + 1)
+                used_prefixes.append(prefix)
+                prefix += "_"
+                # Update the parameters of this constituent model to use the
+                # new prefix
+                for par in submodel_pars:
+                    par.id = prefix + par.id[2:]
+                    par.name = prefix + par.name[2:]
+                    if par.length_control is not None:
+                        par.length_control = prefix + par.length_control[2:]
+                i += npars
+    for part in parts:
         # Parameter prefix per model, A_, B_, ...
         # Note that prefix must also be applied to id and length_control
         # to support vector parameters
+        prefix = chr(ord('A')+k) + '_'
+        scale =  Parameter(prefix+'scale', default=1.0,
+                           description="model intensity for " + part.name)
+        combined_pars.append(scale)
+        prefix = ''
+        if not part.composition:
+            # Model isn't a composition model, so it's parameters don't have a
+            # a prefix. Add the next available prefix
+            prefix = chr(ord('A')+len(used_prefixes))
+            used_prefixes.append(prefix)
+            prefix += '_'
+        if operation == '+':
+            # If model is a sum model, each constituent model gets its own scale parameter
+            scale_prefix = prefix
+            if prefix == '' and part.operation == '*':
+                # `part` is a composition product model. Find the prefixes of
+                # it's parameters to form a new prefix for the scale, eg:
+                # a model with A*B*C will have ABC_scale
+                sub_prefixes = []
+                for param in part.parameters.kernel_parameters:
+                    # Prefix of constituent model
+                    sub_prefix = param.id.split('_')[0]
+                    if sub_prefix not in sub_prefixes:
+                        sub_prefixes.append(sub_prefix)
+                # Concatenate sub_prefixes to form prefix for the scale
+                scale_prefix = ''.join(sub_prefixes) + '_'
+            scale =  Parameter(scale_prefix + 'scale', default=1.0,
+                            description="model intensity for " + part.name)
+            combined_pars.append(scale)
         for p in part.parameters.kernel_parameters:
             p = copy(p)
 …
                 p.length_control = prefix + p.length_control
             combined_pars.append(p)
-        demo.update((prefix+k, v) for k, v in part.demo.items()
-                    if k != "background")
-    #print("pars",combined_pars)
     parameters = ParameterTable(combined_pars)
     parameters.max_pd = sum(part.parameters.max_pd for part in parts)
+    def random():
+        combined_pars = {}
+        for k, part in enumerate(parts):
+            prefix = chr(ord('A')+k) + '_'
+            pars = part.random()
+            combined_pars.update((prefix+k, v) for k, v in pars.items())
+        return combined_pars
     model_info = ModelInfo()
+    model_info.id = '+'.join(part.id for part in parts)
+    model_info.name = ' + '.join(part.name for part in parts)
+    model_info.id = operation.join(part.id for part in parts)
+    model_info.operation = operation
+    model_info.name = '(' + operation.join(part.name for part in parts) + ')'
     model_info.filename = None
     model_info.title = 'Mixture model with ' + model_info.name
 …
     model_info.category = "custom"
     model_info.parameters = parameters
+    model_info.random = random
     #model_info.single = any(part['single'] for part in parts)
     model_info.structure_factor = False
 …
     # Remember the component info blocks so we can build the model
     model_info.composition = ('mixture', parts)
-    model_info.demo = demo
     return model_info
 …
         self.info = model_info
         self.parts = parts
+        self.dtype = parts[0].dtype
     def make_kernel(self, q_vectors):
 …
         self.kernels = kernels
         self.dtype = self.kernels[0].dtype
+        self.operation = model_info.operation
         self.results = []  # type: List[np.ndarray]
 …
         # remember the parts for plotting later
         self.results = []  # type: List[np.ndarray]
-        offset = 2 # skip scale & background
         parts = MixtureParts(self.info, self.kernels, call_details, values)
         for kernel, kernel_details, kernel_values in parts:
             #print("calling kernel", kernel.info.name)
             result = kernel(kernel_details, kernel_values, cutoff, magnetic)
+            #print(kernel.info.name, result)
+            total += result
+            result = np.array(result).astype(kernel.dtype)
+            # print(kernel.info.name, result)
+            if self.operation == '+':
+                total += result
+            elif self.operation == '*':
+                if np.all(total) == 0.0:
+                    total = result
+                else:
+                    total *= result
             self.results.append(result)
 …
         self.part_num += 1
+        self.par_index += info.parameters.npars + 1
+        self.par_index += info.parameters.npars
+        if self.model_info.operation == '+':
+            self.par_index += 1 # Account for each constituent model's scale param
         self.mag_index += 3 * len(info.parameters.magnetism_index)
 …
         # which includes the initial scale and background parameters.
         # We want the index into the weight length/offset for each parameter.
+        # Exclude the initial scale and background, so subtract two, but each
+        # component has its own scale factor which we need to skip when
+        # constructing the details for the kernel, so add one, giving a
+        # net subtract one.
+        index = slice(par_index - 1, par_index - 1 + info.parameters.npars)
+        # Exclude the initial scale and background, so subtract two. If we're
+        # building an addition model, each component has its own scale factor
+        # which we need to skip when constructing the details for the kernel, so
+        # add one, giving a net subtract one.
+        diff = 1 if self.model_info.operation == '+' else 2
+        index = slice(par_index - diff, par_index - diff + info.parameters.npars)
         length = full.length[index]
         offset = full.offset[index]
 …
     def _part_values(self, info, par_index, mag_index):
         # type: (ModelInfo, int, int) -> np.ndarray
+        #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1])
+        scale = self.values[par_index]
+        pars = self.values[par_index + 1:par_index + info.parameters.npars + 1]
+        # Set each constituent model's scale to 1 if this is a multiplication model
+        scale = self.values[par_index] if self.model_info.operation == '+' else 1.0
+        diff = 1 if self.model_info.operation == '+' else 0 # Skip scale if addition model
+        pars = self.values[par_index + diff:par_index + info.parameters.npars + diff]
         nmagnetic = len(info.parameters.magnetism_index)
         if nmagnetic:

sasmodels/model_test.py

-                      rbedb9b0
+                      r65314f7
                 ({}, 'VR', None),
+                ]
+            tests = smoke_tests + self.info.tests
+            tests = smoke_tests
+            if self.info.tests is not None:
+                tests += self.info.tests
             try:
                 model = build_model(self.info, dtype=self.dtype,
 …
         stream.writeln(traceback.format_exc())
         return
     # Run the test suite
     suite.run(result)

sasmodels/modelinfo.py

-                      r724257c
+                      ra85a569
                          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 __getitem__(self, key):
+        # Find the parameter definition
+        for par in self.call_parameters:
+            if par.name == key:
+                break
+        else:
+            raise KeyError("unknown parameter %r"%key)
+        return par
     def _set_vector_lengths(self):
 …
     models when the model is first called, not when the model is loaded.
     """
+    if hasattr(kernel_module, "model_info"):
+        # Custom sum/multi models
+        return kernel_module.model_info
     info = ModelInfo()
     #print("make parameter table", kernel_module.parameters)
 …
     info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
     info.single = getattr(kernel_module, 'single', not callable(info.Iq))
+    info.random = getattr(kernel_module, 'random', None)
     # multiplicity info

sasmodels/models/adsorbed_layer.py

-                      rb0c4271
+                      r8f04da4
 Iq.vectorized = True  # Iq accepts an array of q values
+def random():
+    # only care about the value of second_moment:
+    #    curve = scale * e**(-second_moment^2 q^2)/q^2
+    #    scale = 6 pi/100 (contrast/density*absorbed_amount)^2 * Vf/radius
+    # the remaining parameters can be randomly generated from zero to
+    # twice the default value as done by default in compare.py
+    import numpy as np
+    pars = dict(
+        scale=1,
+        second_moment=10**np.random.uniform(1, 3),
+    )
+    return pars
 # unit test values taken from SasView 3.1.2
 tests = [

sasmodels/models/barbell.py

-                      r9802ab3
+                      r31df0c9
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"]
+def random():
+    import numpy as np
+    # TODO: increase volume range once problem with bell radius is fixed
+    # The issue is that bell radii of more than about 200 fail at high q
+    V = 10**np.random.uniform(7, 9)
+    bar_volume = 10**np.random.uniform(-4, -1)*V
+    bell_volume = V - bar_volume
+    bell_radius = (bell_volume/6)**0.3333  # approximate
+    min_bar = bar_volume/np.pi/bell_radius**2
+    bar_length = 10**np.random.uniform(0, 3)*min_bar
+    bar_radius = np.sqrt(bar_volume/bar_length/np.pi)
+    if bar_radius > bell_radius:
+        bell_radius, bar_radius = bar_radius, bell_radius
+    pars = dict(
+        #background=0,
+        radius_bell=bell_radius,
+        radius=bar_radius,
+        length=bar_length,
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1, background=0,

sasmodels/models/bcc_paracrystal.py

-                      r69e1afc
+                      r8f04da4
 .. figure:: img/parallelepiped_angle_definition.png
     Orientation of the crystal with respect to the scattering plane, when
+    Orientation of the crystal with respect to the scattering plane, when
     $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis).
 …
 source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"]
+# parameters for demo
+demo = dict(
+    scale=1, background=0,
+    dnn=220, d_factor=0.06, sld=4, sld_solvent=1,
+    radius=40,
+    theta=60, phi=60, psi=60,
+    radius_pd=.2, radius_pd_n=2,
+    theta_pd=15, theta_pd_n=0,
+    phi_pd=15, phi_pd_n=0,
+    psi_pd=15, psi_pd_n=0,
+def random():
+    import numpy as np
+    # Define lattice spacing as a multiple of the particle radius
+    # using the formulat a = 4 r/sqrt(3).  Systems which are ordered
+    # are probably mostly filled, so use a distribution which goes from
+    # zero to one, but leaving 90% of them within 80% of the
+    # maximum bcc packing.  Lattice distortion values are empirically
+    # useful between 0.01 and 0.7.  Use an exponential distribution
+    # in this range 'cuz its easy.
+    radius = 10**np.random.uniform(1.3, 4)
+    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
+    dnn_fraction = np.random.beta(a=10, b=1)
+    dnn = radius*4/np.sqrt(3)/dnn_fraction
+    pars = dict(
+        #sld=1, sld_solvent=0, scale=1, background=1e-32,
+        dnn=dnn,
+        d_factor=d_factor,
+        radius=radius,
+    )
+    return pars
 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
 # add 2d test later
 q =4.*pi/220.
+q = 4.*pi/220.
 tests = [
+    [{ },
+     [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]],
+    [{'theta':20.0,'phi':30,'psi':40.0},(-0.017,0.035),2082.20264399 ],
+    [{'theta':20.0,'phi':30,'psi':40.0},(-0.081,0.011),0.436323144781 ]
+    [{}, [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]],
+    [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.017, 0.035), 2082.20264399],
+    [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.081, 0.011), 0.436323144781],
+    ]

sasmodels/models/be_polyelectrolyte.py

-                      r3330bb4
+                      r8f04da4
     r_{0}^2 = \frac{1}{\alpha \sqrt{C_a} \left( b/\sqrt{48\pi L_b}\right)}
 where
+where
 $K$ is the contrast factor for the polymer which is defined differently than in
 …
     a = b_p - (v_p/v_s) b_s
 where $b_p$ and $b_s$ are sum of the scattering lengths of the atoms
 constituting the monomer of the polymer and the sum of the scattering lengths
 …
 $L_b$ is the Bjerrum length(|Ang|) - **Note:** This parameter needs to be
 kept constant for a given solvent and temperature!
+kept constant for a given solvent and temperature!
 $h$ is the virial parameter (|Ang^3|/mol) - **Note:** See [#Borue]_ for the
 …
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
 * **Last Modified by:** Paul Kienzle **Date:** July 24, 2016
 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:**
+* **Last Reviewed by:** Paul Butler and Richard Heenan **Date:**
   October 07, 2016
 """
 …
 Iq.vectorized = True  # Iq accepts an array of q values
+def random():
+    import numpy as np
+    # TODO: review random be_polyelectrolyte model generation
+    pars = dict(
+        scale=10000, #background=0,
+        #polymer_concentration=0.7,
+        polymer_concentration=np.random.beta(5, 3), # around 70%
+        #salt_concentration=0.0,
+        # keep salt concentration extremely low
+        # and use explicit molar to match polymer concentration
+        salt_concentration=np.random.beta(1, 100)*6.022136e-4,
+        #contrast_factor=10.0,
+        contrast_fact=np.random.uniform(1, 100),
+        #bjerrum_length=7.1,
+        bjerrum_length=np.random.uniform(1, 10),
+        #virial_param=12.0,
+        virial_param=np.random.uniform(-1000, 30),
+        #monomer_length=10.0,
+        monomer_length=10.0**(4*np.random.beta(1.5, 3)),
+        #ionization_degree=0.05,
+        ionization_degree=np.random.beta(1.5, 4),
+    )
+    return pars
 demo = dict(scale=1, background=0.1,

sasmodels/models/binary_hard_sphere.py

-                      r925ad6e
+                      r30b60d2
     :nowrap:
     \begin{align}
+    \begin{align*}
     x &= \frac{(\phi_2 / \phi)\alpha^3}{(1-(\phi_2/\phi) + (\phi_2/\phi)
     \alpha^3)} \\
     \phi &= \phi_1 + \phi_2 = \text{total volume fraction} \\
     \alpha &= R_1/R_2 = \text{size ratio}
     \end{align}
+    \end{align*}
 The 2D scattering intensity is the same as 1D, regardless of the orientation of
 …
 source = ["lib/sas_3j1x_x.c", "binary_hard_sphere.c"]
+def random():
+    import numpy as np
+    # TODO: binary_hard_sphere fails at low qr
+    radius_lg = 10**np.random.uniform(2, 4.7)
+    radius_sm = 10**np.random.uniform(2, 4.7)
+    volfraction_lg = 10**np.random.uniform(-3, -0.3)
+    volfraction_sm = 10**np.random.uniform(-3, -0.3)
+    # TODO: Get slightly different results if large and small are swapped
+    # modify the model so it doesn't care which is which
+    if radius_lg < radius_sm:
+        radius_lg, radius_sm = radius_sm, radius_lg
+        volfraction_lg, volfraction_sm = volfraction_sm, volfraction_lg
+    pars = dict(
+        radius_lg=radius_lg,
+        radius_sm=radius_sm,
+        volfraction_lg=volfraction_lg,
+        volfraction_sm=volfraction_sm,
+    )
+    return pars
 # parameters for demo and documentation
 demo = dict(sld_lg=3.5, sld_sm=0.5, sld_solvent=6.36,

sasmodels/models/broad_peak.py

-                      r43fe34b
+                      r0bdddc2
 Iq.vectorized = True  # Iq accepts an array of q values
+def random():
+    import numpy as np
+    pars = dict(
+        scale=1,
+        porod_scale=10**np.random.uniform(-8, -5),
+        porod_exp=np.random.uniform(1, 6),
+        lorentz_scale=10**np.random.uniform(0.3, 6),
+        lorentz_length=10**np.random.uniform(0, 2),
+        peak_pos=10**np.random.uniform(-3, -1),
+        lorentz_exp=np.random.uniform(1, 4),
+    )
+    pars['lorentz_length'] /= pars['peak_pos']
+    pars['lorentz_scale'] *= pars['porod_scale'] / pars['peak_pos']**pars['porod_exp']
+    #pars['porod_scale'] = 0.
+    return pars
 demo = dict(scale=1, background=0,
             porod_scale=1.0e-05, porod_exp=3,

sasmodels/models/capped_cylinder.py

-                      r9802ab3
+                      r31df0c9
 .. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230
 .. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda
+.. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda
    and errata)
 …
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"]
+def random():
+    import numpy as np
+    # TODO: increase volume range once problem with bell radius is fixed
+    # The issue is that bell radii of more than about 200 fail at high q
+    V = 10**np.random.uniform(7, 9)
+    bar_volume = 10**np.random.uniform(-4, -1)*V
+    bell_volume = V - bar_volume
+    bell_radius = (bell_volume/6)**0.3333  # approximate
+    min_bar = bar_volume/np.pi/bell_radius**2
+    bar_length = 10**np.random.uniform(0, 3)*min_bar
+    bar_radius = np.sqrt(bar_volume/bar_length/np.pi)
+    if bar_radius > bell_radius:
+        bell_radius, bar_radius = bar_radius, bell_radius
+    pars = dict(
+        #background=0,
+        radius_cap=bell_radius,
+        radius=bar_radius,
+        length=bar_length,
+    )
+    return pars
 demo = dict(scale=1, background=0,
             sld=6, sld_solvent=1,

sasmodels/models/core_multi_shell.py

-                      r5a0b3d7
+                      r1511c37c
         sld_shell: the SLD of the shell#
         A_shell#: the coefficient in the exponential function
     scale: 1.0 if data is on absolute scale
     volfraction: volume fraction of spheres
 …
 source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"]
+def random():
+    import numpy as np
+    num_shells = np.minimum(np.random.poisson(3)+1, 10)
+    total_radius = 10**np.random.uniform(1.7, 4)
+    thickness = np.random.exponential(size=num_shells+1)
+    thickness *= total_radius/np.sum(thickness)
+    pars = dict(
+        #background=0,
+        n=num_shells,
+        radius=thickness[0],
+    )
+    for k, v in enumerate(thickness[1:]):
+        pars['thickness%d'%(k+1)] = v
+    return pars
 def profile(sld_core, radius, sld_solvent, n, sld, thickness):

sasmodels/models/core_shell_bicelle.py

-                      r9802ab3
+                      r30b60d2
 .. figure:: img/core_shell_bicelle_parameters.png
    Cross section of cylindrical symmetry model used here. Users will have
    to decide how to distribute "heads" and "tails" between the rim, face
+   Cross section of cylindrical symmetry model used here. Users will have
+   to decide how to distribute "heads" and "tails" between the rim, face
    and core regions in order to estimate appropriate starting parameters.
 …
 .. math::
     \rho(r) =
       \begin{cases}
+    \rho(r) =
+      \begin{cases}
       &\rho_c \text{ for } 0 \lt r \lt R; -L \lt z\lt L \\[1.5ex]
       &\rho_f \text{ for } 0 \lt r \lt R; -(L+2t) \lt z\lt -L;
 …
     I(Q,\alpha) = \frac{\text{scale}}{V_t} \cdot
         F(Q,\alpha)^2.sin(\alpha) + \text{background}
+        F(Q,\alpha)^2 \cdot sin(\alpha) + \text{background}
 where
 .. math::
+    :nowrap:
     \begin{align}
     F(Q,\alpha) = &\bigg[
+    \begin{align*}
+    F(Q,\alpha) = &\bigg[
     (\rho_c - \rho_f) V_c \frac{2J_1(QRsin \alpha)}{QRsin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\
     &+(\rho_f - \rho_r) V_{c+f} \frac{2J_1(QRsin\alpha)}{QRsin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} \\
     &+(\rho_r - \rho_s) V_t \frac{2J_1(Q(R+t_r)sin\alpha)}{Q(R+t_r)sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha}
     \bigg]
     \end{align}
+    \end{align*}
 where $V_t$ is the total volume of the bicelle, $V_c$ the volume of the core,
 …
 cylinders is then given by integrating over all possible $\theta$ and $\phi$.
 For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data,
+For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data,
 see the :ref:`cylinder` model for further information.
 Our implementation of the scattering kernel and the 1D scattering intensity
 …
 title = "Circular cylinder with a core-shell scattering length density profile.."
 description = """
     P(q,alpha)= (scale/Vs)*f(q)^(2) + bkg,  where:
+    P(q,alpha)= (scale/Vs)*f(q)^(2) + bkg,  where:
     f(q)= Vt(sld_rim - sld_solvent)* sin[qLt.cos(alpha)/2]
     /[qLt.cos(alpha)/2]*J1(qRout.sin(alpha))
 …
           "core_shell_bicelle.c"]
+def random():
+    import numpy as np
+    pars = dict(
+        radius=10**np.random.uniform(1.3, 3),
+        length=10**np.random.uniform(1.3, 4),
+        thick_rim=10**np.random.uniform(0, 1.7),
+        thick_face=10**np.random.uniform(0, 1.7),
+    )
+    return pars
 demo = dict(scale=1, background=0,
             radius=20.0,

sasmodels/models/core_shell_bicelle_elliptical.py

-                      r9802ab3
+                      r30b60d2
 core-shell scattering length density profile. Thus this is a variation
 of the core-shell bicelle model, but with an elliptical cylinder for the core.
 Outer shells on the rims and flat ends may be of different thicknesses and
+Outer shells on the rims and flat ends may be of different thicknesses and
 scattering length densities. The form factor is normalized by the total particle volume.
 …
 .. figure:: img/core_shell_bicelle_parameters.png
    Cross section of model used here. Users will have
    to decide how to distribute "heads" and "tails" between the rim, face
+   Cross section of model used here. Users will have
+   to decide how to distribute "heads" and "tails" between the rim, face
    and core regions in order to estimate appropriate starting parameters.
 …
 .. math::
     \rho(r) =
       \begin{cases}
+    \rho(r) =
+      \begin{cases}
       &\rho_c \text{ for } 0 \lt r \lt R; -L \lt z\lt L \\[1.5ex]
       &\rho_f \text{ for } 0 \lt r \lt R; -(L+2t) \lt z\lt -L;
 …
     I(Q,\alpha,\psi) = \frac{\text{scale}}{V_t} \cdot
         F(Q,\alpha, \psi)^2.sin(\alpha) + \text{background}
+        F(Q,\alpha, \psi)^2 \cdot sin(\alpha) + \text{background}
 where a numerical integration of $F(Q,\alpha, \psi)^2.sin(\alpha)$ is carried out over \alpha and \psi for:
+where a numerical integration of $F(Q,\alpha, \psi)^2 \cdot sin(\alpha)$ is carried out over \alpha and \psi for:
 .. math::
+    :nowrap:
         \begin{align}
     F(Q,\alpha,\psi) = &\bigg[
+    \begin{align*}
+    F(Q,\alpha,\psi) = &\bigg[
     (\rho_c - \rho_f) V_c \frac{2J_1(QR'sin \alpha)}{QR'sin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\
     &+(\rho_f - \rho_r) V_{c+f} \frac{2J_1(QR'sin\alpha)}{QR'sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} \\
     &+(\rho_r - \rho_s) V_t \frac{2J_1(Q(R'+t_r)sin\alpha)}{Q(R'+t_r)sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha}
     \bigg]
     \end{align}
+    \end{align*}
 where
 …
     R'=\frac{R}{\sqrt{2}}\sqrt{(1+X_{core}^{2}) + (1-X_{core}^{2})cos(\psi)}
 and $V_t = \pi.(R+t_r)(Xcore.R+t_r)^2.(L+2.t_f)$ is the total volume of the bicelle,
 $V_c = \pi.Xcore.R^2.L$ the volume of the core, $V_{c+f} = \pi.Xcore.R^2.(L+2.t_f)$
+and $V_t = \pi.(R+t_r)(Xcore.R+t_r)^2.(L+2.t_f)$ is the total volume of the bicelle,
+$V_c = \pi.Xcore.R^2.L$ the volume of the core, $V_{c+f} = \pi.Xcore.R^2.(L+2.t_f)$
 the volume of the core plus the volume of the faces, $R$ is the radius
 of the core, $Xcore$ is the axial ratio of the core, $L$ the length of the core,
 $t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the usual
 first order bessel function. The core has radii $R$ and $Xcore.R$ so is circular,
 as for the core_shell_bicelle model, for $Xcore$ =1. Note that you may need to
 limit the range of $Xcore$, especially if using the Monte-Carlo algorithm, as
+of the core, $Xcore$ is the axial ratio of the core, $L$ the length of the core,
+$t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the usual
+first order bessel function. The core has radii $R$ and $Xcore.R$ so is circular,
+as for the core_shell_bicelle model, for $Xcore$ =1. Note that you may need to
+limit the range of $Xcore$, especially if using the Monte-Carlo algorithm, as
 setting radius to $R/Xcore$ and axial ratio to $1/Xcore$ gives an equivalent solution!
 …
 bicelles is then given by integrating over all possible $\alpha$ and $\psi$.
 For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
+For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
 see the :ref:`elliptical-cylinder` model for further information.
 …
 .. figure:: img/elliptical_cylinder_angle_definition.png
     Definition of the angles for the oriented core_shell_bicelle_elliptical particles.
+    Definition of the angles for the oriented core_shell_bicelle_elliptical particles.
 …
 description = """
     core_shell_bicelle_elliptical
     Elliptical cylinder core, optional shell on the two flat faces, and shell of
     uniform thickness on its rim (extending around the end faces).
+    Elliptical cylinder core, optional shell on the two flat faces, and shell of
+    uniform thickness on its rim (extending around the end faces).
     Please see full documentation for equations and further details.
     Involves a double numerical integral around the ellipsoid diameter
 …
           "core_shell_bicelle_elliptical.c"]
+demo = dict(scale=1, background=0,
+            radius=30.0,
+            x_core=3.0,
+            thick_rim=8.0,
+            thick_face=14.0,
+            length=50.0,
+            sld_core=4.0,
+            sld_face=7.0,
+            sld_rim=1.0,
+            sld_solvent=6.0,
+            theta=90,
+            phi=0,
+            psi=0)
+def random():
+    import numpy as np
+    outer_major = 10**np.random.uniform(1, 4.7)
+    outer_minor = 10**np.random.uniform(1, 4.7)
+    # Use a distribution with a preference for thin shell or thin core,
+    # limited by the minimum radius. Avoid core,shell radii < 1
+    min_radius = min(outer_major, outer_minor)
+    thick_rim = np.random.beta(0.5, 0.5)*(min_radius-2) + 1
+    radius_major = outer_major - thick_rim
+    radius_minor = outer_minor - thick_rim
+    radius = radius_major
+    x_core = radius_minor/radius_major
+    outer_length = 10**np.random.uniform(1, 4.7)
+    # Caps should be a small percentage of the total length, but at least one
+    # angstrom long.  Since outer length >= 10, the following will suffice
+    thick_face = 10**np.random.uniform(-np.log10(outer_length), -1)*outer_length
+    length = outer_length - thick_face
+    pars = dict(
+        radius=radius,
+        x_core=x_core,
+        thick_rim=thick_rim,
+        thick_face=thick_face,
+        length=length
+    )
+    return pars
 q = 0.1

sasmodels/models/core_shell_cylinder.py

-                      r9b79f29
+                      r9f6823b
     return whole, whole - core
+def random():
+    import numpy as np
+    outer_radius = 10**np.random.uniform(1, 4.7)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    thickness = outer_radius - radius
+    length = np.random.uniform(1, 4.7)
+    pars = dict(
+        radius=radius,
+        thickness=thickness,
+        length=length,
+    )
+    return pars
 demo = dict(scale=1, background=0,
             sld_core=6, sld_shell=8, sld_solvent=1,

sasmodels/models/core_shell_ellipsoid.py

-                      r9802ab3
+                      r30b60d2
 ellipsoid. This may have some undesirable effects if the aspect ratio of the
 ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$
 - which assumes spheres - will not in any case be valid.  Generating a
 custom product model will enable separate effective volume fraction and effective
+- which assumes spheres - will not in any case be valid.  Generating a
+custom product model will enable separate effective volume fraction and effective
 radius in the $S(q)$.
 …
 .. math::
+    \begin{align}
+    :nowrap:
+    \begin{align*}
     F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\
     &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha)
     \end{align}
+    \end{align*}
 where
 .. math::
 …
    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}
 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
+For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,
 see the :ref:`elliptical-cylinder` model for further information.
 …
     return ellipsoid_ER(polar_outer, equat_outer)
+demo = dict(scale=0.05, background=0.001,
+            radius_equat_core=20.0,
+            x_core=3.0,
+            thick_shell=30.0,
+            x_polar_shell=1.0,
+            sld_core=2.0,
+            sld_shell=1.0,
+            sld_solvent=6.3,
+            theta=0,
+            phi=0)
+def random():
+    import numpy as np
+    V = 10**np.random.uniform(5, 12)
+    outer_polar = 10**np.random.uniform(1.3, 4)
+    outer_equatorial = np.sqrt(V/outer_polar) # ignore 4/3 pi
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    thickness_polar = np.random.beta(0.5, 0.5)*(outer_polar-2) + 1
+    thickness_equatorial = np.random.beta(0.5, 0.5)*(outer_equatorial-2) + 1
+    radius_polar = outer_polar - thickness_polar
+    radius_equatorial = outer_equatorial - thickness_equatorial
+    x_core = radius_polar/radius_equatorial
+    x_polar_shell = thickness_polar/thickness_equatorial
+    pars = dict(
+        #background=0, sld=0, sld_solvent=1,
+        radius_equat_core=radius_equatorial,
+        x_core=x_core,
+        thick_shell=thickness_equatorial,
+        x_polar_shell=x_polar_shell,
+    )
+    return pars
 q = 0.1

sasmodels/models/core_shell_parallelepiped.c

-                      r92dfe0c
+                      rc69d6d6
 double form_volume(double length_a, double length_b, double length_c,
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c);
 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
 …
             double thick_rim_c, double theta, double phi, double psi);
 double form_volume(double length_a, double length_b, double length_c,
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c)
+{
     //return length_a * length_b * length_c;
     return length_a * length_b * length_c +
 .0 * thick_rim_a * length_b * length_c +
+    return length_a * length_b * length_c +
+.0 * thick_rim_a * length_b * length_c +
 .0 * thick_rim_b * length_a * length_c +
 .0 * thick_rim_c * length_a * length_b;
 …
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
+    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
     const double mu = 0.5 * q * length_b;
     //calculate volume before rescaling (in original code, but not used)
     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
     //double vol = length_a * length_b * length_c +
     //       2.0 * thick_rim_a * length_b * length_c +
+    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    //double vol = length_a * length_b * length_c +
+    //       2.0 * thick_rim_a * length_b * length_c +
     //       2.0 * thick_rim_b * length_a * length_c +
     //       2.0 * thick_rim_c * length_a * length_b;
     // Scale sides by B
     const double a_scaled = length_a / length_b;
     const double c_scaled = length_c / length_b;
+    // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG)
+    // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c,
+    // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions
+    // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)
+    double ta = (a_scaled + 2.0*thick_rim_a)/length_b;
+    double tb = (a_scaled + 2.0*thick_rim_b)/length_b;
+    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
+    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
+    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
     double Vin = length_a * length_b * length_c;
 …
     double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
+    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
+    double Vot = Vin + V1 + V2 + V3;
     // Scale factors (note that drC is not used later)
 …
     const double drhoA = (arim_sld-solvent_sld);
     const double drhoB = (brim_sld-solvent_sld);
+    //const double drC_Vot = (crim_sld-solvent_sld)*Vot;
+    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
     // Precompute scale factors for combining cross terms from the shape
     const double scale23 = drhoA*V1;
     const double scale14 = drhoB*V2;
+    const double scale12 = drho0*Vin - scale23 - scale14;
+    const double scale24 = drhoC*V3;
+    const double scale11 = drho0*Vin;
+    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
     // outer integral (with gauss points), integration limits = 0, 1
 …
         // inner integral (with gauss points), integration limits = 0, 1
         double inner_total = 0.0;
+        double inner_total_crim = 0.0;
         for(int j=0; j<76; j++) {
             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
 …
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
             const double si2 = sas_sinx_x(mu_proj * cos_uu);
+            const double si2 = sas_sinx_x(mu_proj * cos_uu );
             const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
             const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
 …
             // Expression in libCylinder.c (neither drC nor Vot are used)
             const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
+            const double form_crim = scale11*si1*si2;
+            // To note also that in csparallelepiped.cpp, there is a function called
+            // cspkernel, which seems to make more sense and has the following comment:
+            //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
+            //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
+            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
+            // This is the function called by csparallelepiped_analytical_2D_scaled,
+            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c
             //  correct FF : sum of square of phase factors
             inner_total += Gauss76Wt[j] * form * form;
+            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
+        }
         inner_total *= 0.5;
+        inner_total_crim *= 0.5;
         // now sum up the outer integral
         const double si = sas_sinx_x(mu * c_scaled * sigma);
+        outer_total += Gauss76Wt[i] * inner_total * si * si;
+        const double si_crim = sas_sinx_x(mu * tc * sigma);
+        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
+    }
     outer_total *= 0.5;
 …
     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
+    // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me,
+    // but for the moment I let it like this until understanding better the code.
+    // the scaling by B.
     double ta = length_a + 2.0*thick_rim_a;
     double tb = length_a + 2.0*thick_rim_b;
     double tc = length_a + 2.0*thick_rim_c;
+    double tb = length_b + 2.0*thick_rim_b;
+    double tc = length_c + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
     double siA = sas_sinx_x(0.5*q*length_a*xhat);
 …
     double siBt = sas_sinx_x(0.5*q*tb*yhat);
     double siCt = sas_sinx_x(0.5*q*tc*zhat);
     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
 …
                + drA*(siAt-siA)*siB*siC*V1
                + drB*siA*(siBt-siB)*siC*V2
                + drC*siA*siB*(siCt*siCt-siC)*V3);
+               + drC*siA*siB*(siCt-siC)*V3);
     return 1.0e-4 * f * f;
+}

sasmodels/models/core_shell_parallelepiped.py

-                      r9b79f29
+                      rfa70e04
 Calculates the form factor for a rectangular solid with a core-shell structure.
 The thickness and the scattering length density of the shell or
+The thickness and the scattering length density of the shell or
 "rim" can be different on each (pair) of faces. However at this time
 the 1D calculation does **NOT** actually calculate a c face rim despite the presence of
 …
 **meaning that there are "gaps" at the corners of the solid.**  Again note that
 $t_C = 0$ currently.
+$t_C = 0$ currently.
 The intensity calculated follows the :ref:`parallelepiped` model, with the
 …
 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$
 and length $(C+2t_C)$ values, after appropriately
 sorting the three dimensions to give an oblate or prolate particle, to give an
+sorting the three dimensions to give an oblate or prolate particle, to give an
 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied.
 …
 title = "Rectangular solid with a core-shell structure."
 description = """
      P(q)=
+     P(q)=
 """
 category = "shape:parallelepiped"
 …
 # VR defaults to 1.0
+def random():
+    import numpy as np
+    outer = 10**np.random.uniform(1, 4.7, size=3)
+    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1
+    length = outer - thick
+    pars = dict(
+        length_a=length[0],
+        length_b=length[1],
+        length_c=length[2],
+        thick_rim_a=thick[0],
+        thick_rim_b=thick[1],
+        thick_rim_c=thick[2],
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1, background=0.0,
 …
 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
+tests = [[{}, 0.2, 0.533149288477],
+         [{}, [0.2], [0.533149288477]],
+         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
+         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
+        ]
+del qx, qy  # not necessary to delete, but cleaner
+if 0:  # pak: model rewrite; need to update tests
+    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
+    tests = [[{}, 0.2, 0.533149288477],
+            [{}, [0.2], [0.533149288477]],
+            [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
+            [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
+            ]
+    del qx, qy  # not necessary to delete, but cleaner

sasmodels/models/core_shell_sphere.py

-                      r925ad6e
+                      r9f6823b
 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
+    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]
 …
     return whole, whole - core
+def random():
+    import numpy as np
+    outer_radius = 10**np.random.uniform(1.3, 4.3)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    thickness = outer_radius - radius
+    pars = dict(
+        radius=radius,
+        thickness=thickness,
+    )
+    return pars
 tests = [
     [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],
      # TODO: VR test suppressed until we sort out new product model
      # and determine what to do with volume ratio.
      #[{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704],
+    # TODO: VR test suppressed until we sort out new product model
+    # and determine what to do with volume ratio.
+    #[{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704],
+     # The SasView test result was 0.00169, with a background of 0.001
+     [{'radius': 60.0, 'thickness': 10.0, 'sld_core': 1.0, 'sld_shell':2.0,
+       'sld_solvent':3.0, 'background':0.0},
+.4, 0.000698838],
+    # The SasView test result was 0.00169, with a background of 0.001
+    [{'radius': 60.0, 'thickness': 10.0, 'sld_core': 1.0, 'sld_shell': 2.0,
+      'sld_solvent': 3.0, 'background': 0.0}, 0.4, 0.000698838],
+]

sasmodels/models/cylinder.py

-                      r9802ab3
+                      r31df0c9
 .. figure:: img/cylinder_angle_definition.png
     Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative
     to the beam line coordinates, plus an indication of their orientation distributions
     which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$
+    Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative
+    to the beam line coordinates, plus an indication of their orientation distributions
+    which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$
     in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes.
 …
     Examples for oriented cylinders.
 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
+The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel
+appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel
 to the $Y$ and $X$ axes of the instrument respectively. Some experimentation may be required to understand the 2d patterns fully.
 (Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such
 situations, with very broad distributions, should still be approached with care.)
+(Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such
+situations, with very broad distributions, should still be approached with care.)
 Validation
 …
     return 0.5 * (ddd) ** (1. / 3.)
+def random():
+    import numpy as np
+    V = 10**np.random.uniform(5, 12)
+    length = 10**np.random.uniform(-2, 2)*V**0.333
+    radius = np.sqrt(V/length/np.pi)
+    pars = dict(
+        #scale=1,
+        #background=0,
+        length=length,
+        radius=radius,
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1, background=0,

sasmodels/models/dab.py

-                      r4962519
+                      r404ebbd
     double numerator   = cube(cor_length);
     double denominator = square(1 + square(q*cor_length));
     return numerator / denominator ;
     """
+def random():
+    import numpy as np
+    pars = dict(
+        scale=10**np.random.uniform(1, 4),
+        cor_length=10**np.random.uniform(0.3, 3),
+        #background = 0,
+    )
+    pars['scale'] /= pars['cor_length']**3
+    return pars
 # ER defaults to 1.0

sasmodels/models/ellipsoid.py

-                      r9b79f29
+                      r92708d8
 # Note: model title and parameter table are inserted automatically
 r"""
 The form factor is normalized by the particle volume
+The form factor is normalized by the particle volume
 Definition
 …
 Plenum Press, New York, 1987.
+A. Isihara. J. Chem. Phys. 18(1950) 1446-1449
 Authorship and Verification
 ----------------------------
 …
 * **Last Modified by:** Paul Kienzle **Date:** March 22, 2017
 """
+from __future__ import division
 from numpy import inf, sin, cos, pi
 …
 def ER(radius_polar, radius_equatorial):
     import numpy as np
 # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
     ee = np.empty_like(radius_polar)
     idx = radius_polar > radius_equatorial
 …
     return 0.5 * ddd ** (1.0 / 3.0)
+def random():
+    import numpy as np
+    V = 10**np.random.uniform(5, 12)
+    radius_polar = 10**np.random.uniform(1.3, 4)
+    radius_equatorial = np.sqrt(V/radius_polar) # ignore 4/3 pi
+    pars = dict(
+        #background=0, sld=0, sld_solvent=1,
+        radius_polar=radius_polar,
+        radius_equatorial=radius_equatorial,
+    )
+    return pars
 demo = dict(scale=1, background=0,

sasmodels/models/elliptical_cylinder.py

-                      r9802ab3
+                      rd9ec8f9
     a = qr'\sin(\alpha)
     b = q\frac{L}{2}\cos(\alpha)
     r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)}
 …
 define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$
 (see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle
 $\Psi$ is the rotational angle around its own long_c axis.
+$\Psi$ is the rotational angle around its own long_c axis.
 All angle parameters are valid and given only for 2D calculation; ie, an
 …
     detector plane, with $\Psi$ = 0.
 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
+The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.
 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the
 cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.)
 The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to
 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation
 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)
+appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the
+cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.)
+The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)
 NB: The 2nd virial coefficient of the cylinder is calculated based on the
 …
 * **Author:**
 * **Last Modified by:**
+* **Last Modified by:**
 * **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016
 …
                            + (length + radius) * (length + pi * radius))
     return 0.5 * (ddd) ** (1. / 3.)
+def random():
+    import numpy as np
+    # V = pi * radius_major * radius_minor * length;
+    V = 10**np.random.uniform(3, 9)
+    length = 10**np.random.uniform(1, 3)
+    axis_ratio = 10**np.random.uniform(0, 2)
+    radius_minor = np.sqrt(V/length/axis_ratio)
+    Vf = 10**np.random.uniform(-4, -2)
+    pars = dict(
+        #background=0, sld=0, sld_solvent=1,
+        scale=1e9*Vf/V,
+        length=length,
+        radius_minor=radius_minor,
+        axis_ratio=axis_ratio,
+    )
+    return pars
 q = 0.1
 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct!

sasmodels/models/fcc_paracrystal.py

-                      r69e1afc
+                      r8f04da4
 .. figure:: img/parallelepiped_angle_definition.png
     Orientation of the crystal with respect to the scattering plane, when
+    Orientation of the crystal with respect to the scattering plane, when
     $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis).
 …
 source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"]
+# parameters for demo
+demo = dict(scale=1, background=0,
+            dnn=220, d_factor=0.06, sld=4, sld_solvent=1,
+            radius=40,
+            theta=60, phi=60, psi=60,
+            radius_pd=.2, radius_pd_n=0.2,
+            theta_pd=15, theta_pd_n=0,
+            phi_pd=15, phi_pd_n=0,
+            psi_pd=15, psi_pd_n=0,
+           )
+def random():
+    import numpy as np
+    # copied from bcc_paracrystal
+    radius = 10**np.random.uniform(1.3, 4)
+    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7
+    dnn_fraction = np.random.beta(a=10, b=1)
+    dnn = radius*4/np.sqrt(2)/dnn_fraction
+    pars = dict(
+        #sld=1, sld_solvent=0, scale=1, background=1e-32,
+        dnn=dnn,
+        d_factor=d_factor,
+        radius=radius,
+    )
+    return pars
 # april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!
 q =4.*pi/220.
+q = 4.*pi/220.
 tests = [
+    [{ },
+     [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]],
+     [{}, (-0.047,-0.007), 238.103096286],
+     [{}, (0.053,0.063), 0.863609587796 ],
+    [{}, [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]],
+    [{}, (-0.047, -0.007), 238.103096286],
+    [{}, (0.053, 0.063), 0.863609587796],
+]

sasmodels/models/flexible_cylinder.py

-                      r42356c8
+                      r2573fa1
 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/wrc_cyl.c", "flexible_cylinder.c"]
+demo = dict(scale=1.0, background=0.0001,
+            length=1000.0,
+            kuhn_length=100.0,
+            radius=20.0,
+            sld=1.0,
+            sld_solvent=6.3)
+def random():
+    import numpy as np
+    length = 10**np.random.uniform(2, 6)
+    radius = 10**np.random.uniform(1, 3)
+    kuhn_length = 10**np.random.uniform(-2, 0)*length
+    pars = dict(
+        length=length,
+        radius=radius,
+        kuhn_length=kuhn_length,
+    )
+    return pars
 tests = [
     # Accuracy tests based on content in test/utest_other_models.py
+    # Currently fails in OCL
+    # [{'length':     1000.0,
+    #  'kuhn_length': 100.0,
+    #  'radius':       20.0,
+    #  'sld':           1.0,
+    #  'sld_solvent':   6.3,
+    #  'background':    0.0001,
+    #  }, 0.001, 3509.2187],
+    [{'length':     1000.0,  # test T1
+      'kuhn_length': 100.0,
+      'radius':       20.0,
+      'sld':           1.0,
+      'sld_solvent':   6.3,
+      'background':    0.0001,
+     }, 0.001, 3509.2187],
     # Additional tests with larger range of parameters
     [{'length':    1000.0,
+    [{'length':    1000.0,  # test T2
       'kuhn_length': 100.0,
       'radius':       20.0,
 …
       'background':    0.0001,
      }, 1.0, 0.000595345],
     [{'length':        10.0,
+    [{'length':        10.0,  # test T3
       'kuhn_length': 800.0,
       'radius':        2.0,
 …
       'background':    0.001,
      }, 0.1, 1.55228],
     [{'length':        100.0,
+    [{'length':        100.0,  # test T4
       'kuhn_length': 800.0,
       'radius':       50.0,
 …
+    ]
+# There are a few branches in the code that ought to have test values:
+#
+# For length > 4 * kuhn_length
+#        if length > 10 * kuhn_length then C is scaled by 3.06 (L/b)^(-0.44)
+#        q*kuhn_length <= 3.1  => Sexv_new
+#           dS/dQ < 0 has different behaviour from dS/dQ >= 0
+#  T2    q*kuhn_length > 3.1   => a_long
+#
+# For length <= 4 * kuhn_length
+#        q*kuhn_length <= max(1.9/Rg_short, 3.0)  => Sdebye((q*Rg)^2)
+#           q*Rg < 0.5 uses Pade approx, q*Rg > 1.0 uses math lib
+#  T3,T4 q*kuhn_length > max(1.9/Rg_short, 3.0)   => a_short
+#
+# Note that the transitions between branches may be abrupt.  You can see a
+# several percent change around length=10*kuhn_length and length=4*kuhn_length
+# using the following:
+#
+#    sascomp flexible_cylinder -calc=double -sets=10 length=10*kuhn_length,10.000001*kuhn_length
+#    sascomp flexible_cylinder -calc=double -sets=10 length=4*kuhn_length,4.000001*kuhn_length
+#
+# The transition between low q and high q around q*kuhn_length = 3 seems
+# to be good to 4 digits or better.  This was tested by computing the value
+# on each branches near the transition point and reporting the relative error
+# for kuhn lengths of 10, 100 and 1000 and a variety of length:kuhn_length
+# ratios.

sasmodels/models/flexible_cylinder_elliptical.py

-                      r40a87fa
+                      r31df0c9
           "flexible_cylinder_elliptical.c"]
+demo = dict(scale=1.0, background=0.0001,
+            length=1000.0,
+            kuhn_length=100.0,
+            radius=20.0,
+            axis_ratio=1.5,
+            sld=1.0,
+            sld_solvent=6.3)
+def random():
+    import numpy as np
+    length = 10**np.random.uniform(2, 6)
+    radius = 10**np.random.uniform(1, 3)
+    axis_ratio = 10**np.random.uniform(-1, 1)
+    kuhn_length = 10**np.random.uniform(-2, -0.7)*length  # at least 10 segments
+    pars = dict(
+        length=length,
+        radius=radius,
+        axis_ratio=axis_ratio,
+        kuhn_length=kuhn_length,
+    )
+    return pars
 tests = [

sasmodels/models/fractal.py

-                      rdf89d77
+                      r1511c37c
 where $\xi$ is the correlation length representing the cluster size and $D_f$
 is the fractal dimension, representing the self similarity of the structure.
 Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function
+is the fractal dimension, representing the self similarity of the structure.
+Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function
 diverges at $D_f=0$ and $D_f=1$.
 …
 """
+from __future__ import division
 from numpy import inf
 …
 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"]
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(0.7, 4)
+    #radius = 5
+    cor_length = 10**np.random.uniform(0.7, 2)*radius
+    #cor_length = 20*radius
+    volfraction = 10**np.random.uniform(-3, -1)
+    #volfraction = 0.05
+    fractal_dim = 2*np.random.beta(3, 4) + 1
+    #fractal_dim = 2
+    pars = dict(
+        #background=0, sld_block=1, sld_solvent=0,
+        volfraction=volfraction,
+        radius=radius,
+        cor_length=cor_length,
+        fractal_dim=fractal_dim,
+    )
+    return pars
 demo = dict(volfraction=0.05,
             radius=5.0,

sasmodels/models/fractal_core_shell.py

-                      r925ad6e
+                      rca04add
     \frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3}+
 V_s(\rho_s-\rho_{solv})
+    \frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right]^2
+    S(q) &= 1 + \frac{D_f\ \Gamma\!(D_f-1)}{[1+1/(q\xi)^2]^{(D_f-1)/2}}
+    \frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right]^2 \\
+    S(q) &= 1 + \frac{D_f\ \Gamma\!(D_f-1)}{[1+1/(q\xi)^2]^{(D_f-1)/2}}
     \frac{\sin[(D_f-1)\tan^{-1}(q\xi)]}{(qr_s)^{D_f}}
 …
 $\rho_{solv}$ are the scattering length densities of the core, shell, and
 solvent respectively, $r_c$ and $r_s$ are the radius of the core and the radius
 of the whole particle respectively, $D_f$ is the fractal dimension, and |xi| the
+of the whole particle respectively, $D_f$ is the fractal dimension, and $\xi$ the
 correlation length.
 Polydispersity of radius and thickness are also provided for.
 …
           "lib/fractal_sq.c", "fractal_core_shell.c"]
+def random():
+    import numpy as np
+    outer_radius = 10**np.random.uniform(0.7, 4)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    thickness = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    radius = outer_radius - thickness
+    cor_length = 10**np.random.uniform(0.7, 2)*outer_radius
+    volfraction = 10**np.random.uniform(-3, -1)
+    fractal_dim = 2*np.random.beta(3, 4) + 1
+    pars = dict(
+        #background=0, sld_block=1, sld_solvent=0,
+        volfraction=volfraction,
+        radius=radius,
+        cor_length=cor_length,
+        fractal_dim=fractal_dim,
+    )
+    return pars
 demo = dict(scale=0.05,
             background=0,

sasmodels/models/fuzzy_sphere.py

-                      r925ad6e
+                      r31df0c9
 # VR defaults to 1.0
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(1, 4.7)
+    fuzziness = 10**np.random.uniform(-2, -0.5)*radius  # 1% to 31% fuzziness
+    pars = dict(
+        radius=radius,
+        fuzziness=fuzziness,
+    )
+    return pars
 demo = dict(scale=1, background=0.001,
             sld=1, sld_solvent=3,

sasmodels/models/gauss_lorentz_gel.py

-                      ra807206
+                      r48462b0
 but typically a physical rather than chemical network.
 It is modeled as a sum of a low-q exponential decay (which happens to
 give a functional form similar to Guinier scattering, so interpret with
+give a functional form similar to Guinier scattering, so interpret with
 care) plus a Lorentzian at higher-q values. See also the gel_fit model.
 …
+def random():
+    import numpy as np
+    gauss_scale = 10**np.random.uniform(1, 3)
+    lorentz_scale = 10**np.random.uniform(1, 3)
+    cor_length_static = 10**np.random.uniform(0, 3)
+    cor_length_dynamic = 10**np.random.uniform(0, 3)
+    pars = dict(
+        #background=0,
+        scale=1,
+        gauss_scale=gauss_scale,
+        lorentz_scale=lorentz_scale,
+        cor_length_static=cor_length_static,
+        cor_length_dynamic=cor_length_dynamic,
+    )
+    return pars
 demo = dict(scale=1, background=0.1,
             gauss_scale=100.0,

sasmodels/models/gaussian_peak.py

-                      ra807206
+                      r48462b0
 # VR defaults to 1.0
+demo = dict(scale=1, background=0, peak_pos=0.05, sigma=0.005)
+def random():
+    import numpy as np
+    peak_pos = 10**np.random.uniform(-3, -1)
+    sigma = 10**np.random.uniform(-1.3, -0.3)*peak_pos
+    scale = 10**np.random.uniform(0, 4)
+    pars = dict(
+        #background=1e-8,
+        scale=scale,
+        peak_pos=peak_pos,
+        sigam=sigma,
+    )
+    return pars

sasmodels/models/gel_fit.c

-                      ra807206
+                      r48462b0
     double lorentzian_term = square(q*cor_length);
     lorentzian_term = 1.0 + ((fractal_dim + 1.0)/3.0)*lorentzian_term;
     lorentzian_term = pow(lorentzian_term, (fractal_dim/2.0) );
+    lorentzian_term = pow(lorentzian_term, fractal_dim/2.0 );
     // Exponential Term
     ////////////////////////double d(x[i]*x[i]*rg*rg);
     double exp_term = square(q*rg);
     exp_term = exp(-1.0*(exp_term/3.0));
+    exp_term = exp(-exp_term/3.0);
     // Scattering Law

sasmodels/models/gel_fit.py

-                      ra807206
+                      r48462b0
 source = ["gel_fit.c"]
+def random():
+    import numpy as np
+    guinier_scale = 10**np.random.uniform(1, 3)
+    lorentz_scale = 10**np.random.uniform(1, 3)
+    rg = 10**np.random.uniform(1, 5)
+    fractal_dim = np.random.uniform(0, 6)
+    cor_length = 10**np.random.uniform(0, 3)
+    pars = dict(
+        #background=0,
+        scale=1,
+        guinier_scale=guinier_scale,
+        lorentz_scale=lorentz_scale,
+        rg=rg,
+        fractal_dim=fractal_dim,
+        cor_length=cor_length
+    )
+    return pars
 demo = dict(background=0.01,
             guinier_scale=1.7,

sasmodels/models/guinier.py

-                      r3330bb4
+                      r48462b0
 description = """
  I(q) = scale.exp ( - rg^2 q^2 / 3.0 )
     List of default parameters:
     scale = scale
 …
 """
+def random():
+    import numpy as np
+    scale = 10**np.random.uniform(1, 4)
+    # Note: compare.py has Rg cutoff of 1e-30 at q=1 for guinier, so use that
+    # log_10 Ae^(-(q Rg)^2/3) = log_10(A) - (q Rg)^2/ (3 ln 10) > -30
+    #   => log_10(A) > Rg^2/(3 ln 10) - 30
+    q_max = 1.0
+    rg_max = np.sqrt(90*np.log(10) + 3*np.log(scale))/q_max
+    rg = 10**np.random.uniform(0, np.log10(rg_max))
+    pars = dict(
+        #background=0,
+        scale=scale,
+        rg=rg,
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1.0, rg=60.0)

sasmodels/models/guinier_porod.py

-                      rcdcebf1
+                      r48462b0
 and dimensionality of scattering objects, including asymmetric objects
 such as rods or platelets, and shapes intermediate between spheres
 and rods or between rods and platelets, and overcomes some of the
+and rods or between rods and platelets, and overcomes some of the
 deficiencies of the (Beaucage) Unified_Power_Rg model (see Hammouda, 2010).
 …
                         scale = Guinier Scale
                         s = Dimension Variable
                         Rg = Radius of Gyration [A]
+                        Rg = Radius of Gyration [A]
                         porod_exp = Porod Exponent
                         background  = Background [1/cm]"""
 …
 Iq.vectorized = True # Iq accepts an array of q values
+def random():
+    import numpy as np
+    rg = 10**np.random.uniform(1, 5)
+    s = np.random.uniform(0, 3)
+    porod_exp = s + np.random.uniform(0, 3)
+    pars = dict(
+        #scale=1, background=0,
+        rg=rg,
+        s=s,
+        porod_exp=porod_exp,
+    )
+    return pars
 demo = dict(scale=1.5, background=0.5, rg=60, s=1.0, porod_exp=3.0)

sasmodels/models/hardsphere.py

-                      r40a87fa
+                      r8f04da4
 Iq = r"""
       double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH;
       // these are c compiler instructions, can also put normal code inside the "if else" structure
+      // these are c compiler instructions, can also put normal code inside the "if else" structure
       #if FLOAT_SIZE > 4
       // double precision    orig had 0.2, don't call the variable cutoff as PAK already has one called that! Must use UPPERCASE name please.
 …
       #else
       // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good
       #define CUTOFFHS 0.4
+      #define CUTOFFHS 0.4
       #endif
 …
       if(X < CUTOFFHS) {
       // RKH Feb 2016, use Taylor series expansion for small X
       // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures.
       // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient
       // for 5 or 6 significant figures, but I put the X^4 one in anyway
+      // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures.
+      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient
+      // for 5 or 6 significant figures, but I put the X^4 one in anyway
             //FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4;
             // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec)
 …
    """
+def random():
+    import numpy as np
+    pars = dict(
+        scale=1, background=0,
+        radius_effective=10**np.random.uniform(1, 4),
+        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
+    )
+    return pars
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/hayter_msa.py

-                      r5702f52
+                      r8f04da4
         Routine takes absolute value of charge, use HardSphere if charge
         goes to zero.
         In sasview the effective radius and volume fraction may be calculated
+        In sasview the effective radius and volume fraction may be calculated
         from the parameters used in P(Q).
 """
 …
 # ER defaults to 0.0
 # VR defaults to 1.0
+def random():
+    import numpy as np
+    # TODO: too many failures for random hayter_msa parameters
+    pars = dict(
+        scale=1, background=0,
+        radius_effective=10**np.random.uniform(1, 4.7),
+        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction
+        charge=min(int(10**np.random.uniform(0, 1.3)+0.5), 200),
+        temperature=10**np.random.uniform(0, np.log10(450)), # max T = 450
+        #concentration_salt=10**np.random.uniform(-3, 1),
+        dialectconst=10**np.random.uniform(0, 6),
+        #charge=10,
+        #temperature=318.16,
+        concentration_salt=0.0,
+        #dielectconst=71.08,
+    )
+    return pars
 # default parameter set,  use  compare.sh -midQ -linear

sasmodels/models/hollow_cylinder.py

-                      rf102a96
+                      r8f04da4
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
 * **Last Modified by:** Richard Heenan **Date:** October 06, 2016
+* **Last Modified by:** Richard Heenan **Date:** October 06, 2016
    (reparametrised to use thickness, not outer radius)
 * **Last Reviewed by:** Richard Heenan **Date:** October 06, 2016
 …
     return vol_shell, vol_total
+def random():
+    import numpy as np
+    length = 10**np.random.uniform(1, 4.7)
+    outer_radius = 10**np.random.uniform(1, 4.7)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    thickness = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    radius = outer_radius - thickness
+    pars = dict(
+        length=length,
+        radius=radius,
+        thickness=thickness,
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1.0, background=0.0, length=400.0, radius=20.0,

sasmodels/models/hollow_rectangular_prism.py

-                      rab2aea8
+                      r30b60d2
   :nowrap:
   \begin{align}
+  \begin{align*}
   A_{P\Delta}(q) & =  A B C
     \left[\frac{\sin \bigl( q \frac{C}{2} \cos\theta \bigr)}
 …
     \left[ \frac{\sin \bigl[ q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi \bigr]}
     {q \bigl(\frac{B}{2}-\Delta\bigr) \sin\theta \cos\phi} \right]
   \end{align}
+  \end{align*}
 where $A$, $B$ and $C$ are the external sides of the parallelepiped fulfilling
 …
+def random():
+    import numpy as np
+    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
+    # Thickness is limited to 1/2 the smallest dimension
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    min_dim = 0.5*min(a, b, c)
+    thickness = np.random.beta(0.5, 0.5)*(min_dim-2) + 1
+    #print(a, b, c, thickness, thickness/min_dim)
+    pars = dict(
+        length_a=a,
+        b2a_ratio=b/a,
+        c2a_ratio=c/a,
+        thickness=thickness,
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1, background=0,

sasmodels/models/hollow_rectangular_prism_thin_walls.py

-                      rab2aea8
+                      r31df0c9
+def random():
+    import numpy as np
+    a, b, c = 10**np.random.uniform(1, 4.7, size=3)
+    pars = dict(
+        length_a=a,
+        b2a_ratio=b/a,
+        c2a_ratio=c/a,
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1, background=0,

sasmodels/models/lamellar.py

-                      r40a87fa
+                      r1511c37c
     """
+def random():
+    import numpy as np
+    thickness = 10**np.random.uniform(1, 4)
+    pars = dict(
+        thickness=thickness,
+    )
+    return pars
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lamellar_hg.py

-                      ra807206
+                      r1511c37c
     """
+def random():
+    import numpy as np
+    thickness = 10**np.random.uniform(1, 4)
+    length_head = thickness * np.random.uniform(0, 1)
+    length_tail = thickness - length_head
+    pars = dict(
+        length_head=length_head,
+        length_tail=length_tail,
+    )
+    return pars
 # ER defaults to 0.0
 # VR defaults to 1.0

sasmodels/models/lamellar_hg_stack_caille.py

-                      ra57b31d
+                      r1511c37c
 # VR defaults to 1.0
+def random():
+    import numpy as np
+    total_thickness = 10**np.random.uniform(2, 4.7)
+    Nlayers = np.random.randint(2, 200)
+    d_spacing = total_thickness / Nlayers
+    thickness = d_spacing * np.random.uniform(0, 1)
+    length_head = thickness * np.random.uniform(0, 1)
+    length_tail = thickness - length_head
+    Caille_parameter = np.random.uniform(0, 0.8)
+    pars = dict(
+        length_head=length_head,
+        length_tail=length_tail,
+        Nlayers=Nlayers,
+        d_spacing=d_spacing,
+        Caille_parameter=Caille_parameter,
+    )
+    return pars
 demo = dict(
     scale=1, background=0,

sasmodels/models/lamellar_stack_caille.py

-                      ra57b31d
+                      r1511c37c
 source = ["lamellar_stack_caille.c"]
+def random():
+    import numpy as np
+    total_thickness = 10**np.random.uniform(2, 4.7)
+    Nlayers = np.random.randint(2, 200)
+    d_spacing = total_thickness / Nlayers
+    thickness = d_spacing * np.random.uniform(0, 1)
+    Caille_parameter = np.random.uniform(0, 0.8)
+    pars = dict(
+        thickness=thickness,
+        Nlayers=Nlayers,
+        d_spacing=d_spacing,
+        Caille_parameter=Caille_parameter,
+    )
+    return pars
 # No volume normalization despite having a volume parameter
 # This should perhaps be volume normalized?

sasmodels/models/lamellar_stack_paracrystal.py

-                      ra0168e8
+                      r8f04da4
 form_volume = """
     return 1.0;
+    """
+"""
+def random():
+    import numpy as np
+    total_thickness = 10**np.random.uniform(2, 4.7)
+    Nlayers = np.random.randint(2, 200)
+    d_spacing = total_thickness / Nlayers
+    thickness = d_spacing * np.random.uniform(0, 1)
+    # Let polydispersity peak around 15%; 95% < 0.4; max=100%
+    sigma_d = np.random.beta(1.5, 7)
+    pars = dict(
+        thickness=thickness,
+        Nlayers=Nlayers,
+        d_spacing=d_spacing,
+        sigma_d=sigma_d,
+    )
+    return pars
 # ER defaults to 0.0

sasmodels/models/lib/wrc_cyl.c

-                      r92ce163
+                      r18a2bfc
     Functions for WRC implementation of flexible cylinders
 */
+double Sk_WR(double q, double L, double b);
+static double
+AlphaSquare(double x)
+{
+    // Potentially faster. Needs proper benchmarking.
+    // add native_powr to kernel_template
+    //double t = native_powr( (1.0 + (x/3.12)*(x/3.12) +
+    //     (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
+    //return t;
+    return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0);
+}
+//
+static double
+Rgsquarezero(double q, double L, double b)
+{
+    const double r = b/L;
+    return (L*b/6.0) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r))));
+}
+//
+static double
+Rgsquareshort(double q, double L, double b)
+{
+    return AlphaSquare(L/b) * Rgsquarezero(q,L,b);
+}
+//
+static double
+Rgsquare(double q, double L, double b)
+{
+    return AlphaSquare(L/b)*L*b/6.0;
+}
+static double
+sech_WR(double x)
+{
+    return(1/cosh(x));
+}
+static double
+a1long(double q, double L, double b, double p1, double p2, double q0)
+{
+    double C;
+    if( L/b > 10.0) {
+        C = 3.06/pow((L/b),0.44);
+static double
+Rgsquare(double L, double b)
+{
+    const double x = L/b;
+    // Use Horner's method to evaluate Pedersen eq 15:
+    //     alpha^2 = [1.0 + (x/3.12)^2 + (x/8.67)^3] ^ (0.176/3)
+    const double alphasq =
+        pow(1.0 + x*x*(1.027284681130835e-01 + 1.534414548417740e-03*x),
+.866666666666667e-02);
+    return alphasq*L*b/6.0;
+}
+static double
+Rgsquareshort(double L, double b)
+{
+    const double r = b/L;  // = 1/n_b in Pedersen ref.
+    return Rgsquare(L, b) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r))));
+}
+static double
+w_WR(double x)
+{
+    // Pedersen eq. 16:
+    //    w = [1 + tanh((x-C4)/C5)]/2
+    const double C4 = 1.523;
+    const double C5 = 0.1477;
+    return 0.5 + 0.5*tanh((x - C4)/C5);
+}
+static double
+Sdebye(double qsq)
+{
+#if FLOAT_SIZE>4
+#  define DEBYE_CUTOFF 0.25  // 1e-15 error
+#else
+#  define DEBYE_CUTOFF 1.0  // 4e-7 error
+#endif
+    if (qsq < DEBYE_CUTOFF) {
+        const double x = qsq;
+        // mathematica: 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.;
+        const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.;
+        return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.)
+             / ((((B4*x + B3)*x + B2)*x + B1)*x + 1.);
     } else {
         C = 1.0;
+        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq);
+    }
+}
+static double
+a_long(double q, double L, double b/*, double p1, double p2, double q0*/)
+{
+    const double p1 = 4.12;
+    const double p2 = 4.42;
+    const double q0 = 3.1;
+    // Constants C1, ..., C5 derived from least squares fit (Pedersen, eq 13,16)
     const double C1 = 1.22;
     const double C2 = 0.4288;
 …
     const double miu = 0.585;
+    const double Rg2 = Rgsquare(q,L,b);
+    const double Rg22 = Rg2*Rg2;
+    const double Rg = sqrt(Rg2);
+    const double Rgb = Rg*q0/b;
+    const double b2 = b*b;
+    const double C = (L/b>10.0 ? 3.06*pow(L/b, -0.44) : 1.0);
+    //printf("branch B-%d q=%g L=%g b=%g\n", C==1.0, q, L, b);
+    const double r2 = Rgsquare(L,b);
+    const double r = sqrt(r2);
+    const double qr_b = q0*r/b;
+    const double qr_b_sq = qr_b*qr_b;
+    const double qr_b_4 = qr_b_sq*qr_b_sq;
+    const double qr_b_miu = pow(qr_b, -1.0/miu);
+    const double em1_qr_b_sq = expm1(-qr_b_sq);
+    const double sech2 = 1.0/square(cosh((qr_b-C4)/C5));
+    const double w = w_WR(qr_b);
+    const double t1 = pow(q0, 1.0 + p1 + p2)/(b*(p1-p2));
+    const double t2 = C/(15.0*L) * (
+        + 14.0*b*b*em1_qr_b_sq/(q0*qr_b_sq)
+        + 2.0*q0*r2*exp(-qr_b_sq)*(11.0 + 7.0/qr_b_sq));
+    const double t11 = ((C3*qr_b_miu + C2)*qr_b_miu + C1)*qr_b_miu;
+    const double t3 = r*sech2/(2.*C5)*t11;
+    const double t4 = r*(em1_qr_b_sq + qr_b_sq)*sech2 / (C5*qr_b_4);
+    const double t5 = -4.0*r*qr_b*em1_qr_b_sq/qr_b_4 * (1.0 - w);
+    const double t10 = 2.0*(em1_qr_b_sq + qr_b_sq)/qr_b_4 * (1.0 - w); //=Sdebye*(1-w)
+    const double t6 = 4.0*b/q0 * t10;
+    const double t7 = r*((-3.0*C3*qr_b_miu -2.0*C2)*qr_b_miu -1.0*C1)*qr_b_miu/(miu*qr_b);
+    const double t9 = C*b/L * (4.0 - exp(-qr_b_sq) * (11.0 + 7.0/qr_b_sq) + 7.0/qr_b_sq)/15.0;
+    const double t12 = b*b*M_PI/(L*q0*q0) + t2 + t3 - t4 + t5 - t6 + t7*w;
+    const double t13 = -b*M_PI/(L*q0) + t9 + t10 + t11*w;
+    const double a1 = pow(q0,p1)*t13 - t1*pow(q0,-p2)*(t12 + b*p1/q0*t13);
+    const double a2 = t1*pow(q0,-p1)*(t12 + b*p1/q0*t13);
+    const double ans = a1*pow(q*b, -p1) + a2*pow(q*b, -p2) + M_PI/(q*L);
+    return ans;
+}
+static double
+_short(double r2, double exp_qr_b, double L, double b, double p1short,
+        double p2short, double q0)
+{
+    const double qr2 = q0*q0 * r2;
     const double b3 = b*b*b;
+    const double b4 = b3*b;
+    const double q02 = q0*q0;
+    const double q03 = q0*q0*q0;
+    const double q04 = q03*q0;
+    const double q05 = q04*q0;
+    const double Rg02 = Rg2*q02;
+    const double t1 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2))) *
+         ((11.0/15.0 + (7.0*b2)/(15.0*Rg02))) +
+         (7.0*b2)/(15.0*Rg02))));
+    const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*((1.0 + 0.5*(((-1.0) -
+         tanh((-C4 + Rgb/C5)))))));
+    const double t3 = ((C3*pow(Rgb,((-3.0)/miu)) +
+         C2*pow(Rgb,((-2.0)/miu)) +
+         C1*pow(Rgb,((-1.0)/miu))));
+    const double t4 = ((1.0 + tanh(((-C4) + Rgb)/C5)));
+    const double t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) -
+         b*p2*pow(q0,((-1.0) - p1 - p2))));
+    const double t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) +
+         (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) +
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 +
+         (7.0*b2)/(15.0*Rg02)))*Rg2)/b)));
+    const double t7 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) +
+         C2*pow(((Rgb)),((-2.0)/miu)) +
+         C1*pow(((Rgb)),((-1.0)/miu))))*pow(sech_WR(((-C4) +
+         Rgb)/C5),2));
+    const double t8 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*pow(sech_WR(((-C4) + Rgb)/C5),2));
+    const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b -
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-1.0) -
+         tanh(((-C4) + Rgb)/C5))))));
+    const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) +
+         Rgb)/C5))))));
+    const double t11 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) -
+.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),((-1.0) -
+.0/miu)))/miu - (C1*Rg*pow(((Rgb)),((-1.0) -
+.0/miu)))/miu));
+    const double t12 = ((1.0 + tanh(((-C4) + Rgb)/C5)));
+    const double t13 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 +
+          (7.0*b2)/(15.0*q02* Rg2))) +
+          (7.0*b2)/(15.0*Rg02))));
+    const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+          Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) +
+          Rgb)/C5))))));
+    const double t15 = ((C3*pow(((Rgb)),((-3.0)/miu)) +
+        C2*pow(((Rgb)),((-2.0)/miu)) +
+        C1*pow(((Rgb)),((-1.0)/miu))));
+    double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) +
+.5*t3*t4)) + (t5*((pow(q0,(p1 - p2))*
+        (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) -
+        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*t11*t12)) -
+        b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L +
+        t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) +
+        Rgb)/C5)))))))))));
+    return (yy);
+}
+static double
+a2long(double q, double L, double b, double p1, double p2, double q0)
+{
+    double C;
+    if( L/b > 10.0) {
+        C = 3.06/pow((L/b),0.44);
+    } else {
+        C = 1.0;
+    }
+    const double C1 = 1.22;
+    const double C2 = 0.4288;
+    const double C3 = -1.651;
+    const double C4 = 1.523;
+    const double C5 = 0.1477;
+    const double miu = 0.585;
+    const double Rg2 = Rgsquare(q,L,b);
+    const double Rg22 = Rg2*Rg2;
+    const double b2 = b*b;
+    const double b3 = b*b*b;
+    const double b4 = b3*b;
+    const double q02 = q0*q0;
+    const double q03 = q0*q0*q0;
+    const double q04 = q03*q0;
+    const double q05 = q04*q0;
+    const double Rg = sqrt(Rg2);
+    const double Rgb = Rg*q0/b;
+    const double Rg02 = Rg2*q02;
+    const double t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) -
+         b*p2*pow(q0,((-1.0) - p1 - p2)) ));
+    const double t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) +
+         (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) +
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 +
+         (7*b2)/(15.0*Rg02)))*Rg2)/b)))/L;
+    const double t3 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) +
+         C2*pow(((Rgb)),((-2.0)/miu)) +
+         C1*pow(((Rgb)),((-1.0)/miu))))*
+         pow(sech_WR(((-C4) +Rgb)/C5),2.0))/(2.0*C5);
+    const double t4 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*pow(sech_WR(((-C4) +
+         Rgb)/C5),2))/(C5*q04*Rg22);
+    const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b -
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*
+         ((1.0 + 0.5*(((-1.0) - tanh(((-C4) +
+         Rgb)/C5))))))/(q04*Rg22);
+    const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) +
+         Rgb)/C5))))))/(q05*Rg22);
+    const double t7 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) -
+.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),
+         ((-1.0) - 2.0/miu)))/miu - (C1*Rg*pow(((Rgb)),
+         ((-1.0) - 1.0/miu)))/miu));
+    const double t8 = ((1.0 + tanh(((-C4) + Rgb)/C5)));
+    const double t9 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 +
+         (7.0*b2)/(15*Rg02))) + (7.0*b2)/(15.0*Rg02))))/L;
+    const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) +
+          Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) +
+          Rgb)/C5))))))/(q04*Rg22);
+    const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) +
+         t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*
+         (((-((b*M_PI)/(L*q0))) + t9 + t10 +
+.5*((C3*pow(((Rgb)),((-3.0)/miu)) +
+         C2*pow(((Rgb)),((-2.0)/miu)) +
+         C1*pow(((Rgb)),((-1.0)/miu))))*
+         ((1.0 + tanh(((-C4) + Rgb)/C5))))))))));
+    return (yy);
+}
+//
+static double
+a_short(double q, double L, double b, double p1short,
+        double p2short, double factor, double pdiff, double q0)
+{
+    const double Rg2_sh = Rgsquareshort(q,L,b);
+    const double Rg2_sh2 = Rg2_sh*Rg2_sh;
+    const double b3 = b*b*b;
+    const double t1 = ((q0*q0*Rg2_sh)/(b*b));
+    const double Et1 = exp(t1);
+    const double Emt1 = 1.0/Et1;
+    const double q02 = q0*q0;
+    const double q0p = pow(q0,(-4.0 + p1short) );
+    double yy = ((factor/(L*(pdiff)*Rg2_sh2)*
+        ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short +
+.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh -
+.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*M_PI*q02*q0*Rg2_sh2 +
+        Et1*p2short*M_PI*q02*q0*Rg2_sh2))))));
+    return(yy);
+}
+static double
+a1short(double q, double L, double b, double p1short, double p2short, double q0)
+{
+    double factor = 1.0;
+    return a_short(q, L, b, p1short, p2short, factor, p1short - p2short, q0);
+}
+static double
+a2short(double q, double L, double b, double p1short, double p2short, double q0)
+{
+    double factor = -1.0;
+    return a_short(q, L, b, p2short, p1short, factor, p1short-p2short, q0);
+}
+//WR named this w (too generic)
+static double
+w_WR(double x)
+{
+    return 0.5*(1 + tanh((x - 1.523)/0.1477));
+}
+//
+static double
+u1(double q, double L, double b)
+{
+    return Rgsquareshort(q,L,b)*q*q;
+}
+static double
+u_WR(double q, double L, double b)
+{
+    return Rgsquare(q,L,b)*q*q;
+}
+static double
+Sdebye_kernel(double arg)
+{
+    // ORIGINAL
+    double result = 2.0*(exp(-arg) + arg -1.0)/(arg*arg);
+    // CONVERSION 1 from http://herbie.uwplse.org/
+    //
+    // exhibits discontinuity - needs more investigation
+    //double a1 = 1.0/6.0;
+    //double a2 = 1.0/72.0;
+    //double a3 = 1.0/24.0;
+    //double result = pow((1.0 - a1*arg - (a2+a3)*arg*arg), 2);
+    return result;
+}
+static double
+Sdebye(double q, double L, double b)
+{
+    double arg = u_WR(q,L,b);
+    return Sdebye_kernel(arg);
+}
+//
+static double
+Sdebye1(double q, double L, double b)
+{
+    double arg = u1(q,L,b);
+    return Sdebye_kernel(arg);
+}
+//
+    const double q0p = pow(q0, -4.0 + p1short);
+    double yy = 1.0/(L*r2*r2) * b/exp_qr_b*q0p
+        * (8.0*b3*L
+           - 8.0*b3*exp_qr_b*L
+           + 2.0*b3*exp_qr_b*L*p2short
+           - 2.0*b*exp_qr_b*L*p2short*qr2
+           + 4.0*b*exp_qr_b*L*qr2
+           - 2.0*b3*L*p2short
+           + 4.0*b*L*qr2
+           - M_PI*exp_qr_b*qr2*q0*r2
+           + M_PI*exp_qr_b*p2short*qr2*q0*r2);
+    return yy;
+}
+static double
+a_short(double qp, double L, double b
+        /*double p1short, double p2short*/, double q0)
+{
+    const double p1short = 5.36;
+    const double p2short = 5.62;
+    const double r2 = Rgsquareshort(L,b);
+    const double exp_qr_b = exp(r2*square(q0/b));
+    const double pdiff = p1short - p2short;
+    const double a1 = _short(r2,exp_qr_b,L,b,p1short,p2short,q0)/pdiff;
+    const double a2= -_short(r2,exp_qr_b,L,b,p2short,p1short,q0)/pdiff;
+    const double ans = a1*pow(qp*b, -p1short) + a2*pow(qp*b, -p2short) + M_PI/(qp*L);
+    return ans;
+}
 static double
 Sexv(double q, double L, double b)
+{
+    // Pedersen eq 13, corrected by Chen eq A.5, swapping w and 1-w
     const double C1=1.22;
     const double C2=0.4288;
     const double C3=-1.651;
     const double miu = 0.585;
+    const double qRg = q*sqrt(Rgsquare(q,L,b));
+    const double x = pow(qRg, -1.0/miu);
+    double yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) +
+            w_WR(qRg)*x*(C1 + x*(C2 + x*C3));
+    return (yy);
+}
+static double
+Sexvnew(double q, double L, double b)
+{
+    double yy;
+    const double C1 =1.22;
+    const double C2 =0.4288;
+    const double C3 =-1.651;
+    const double miu = 0.585;
+    const double qr = q*sqrt(Rgsquare(L,b));
+    const double qr_miu = pow(qr, -1.0/miu);
+    const double w = w_WR(qr);
+    const double t10 = Sdebye(qr*qr)*(1.0 - w);
+    const double t11 = ((C3*qr_miu + C2)*qr_miu + C1)*qr_miu;
+    return t10 + w*t11;
+}
+// Modified by Yun on Oct. 15,
+static double
+Sexv_new(double q, double L, double b)
+{
+    const double qr = q*sqrt(Rgsquare(L,b));
+    const double qr2 = qr*qr;
+    const double C = (L/b > 10.0) ? 3.06*pow(L/b, -0.44) : 1.0;
+    const double t9 = C*b/L * (4.0 - exp(-qr2) * (11.0 + 7.0/qr2) + 7.0/qr2)/15.0;
+    const double Sexv_orig = Sexv(q, L, b);
+    // calculating the derivative to decide on the correction (cutoff) term?
+    // Note: this is modified from WRs original code
     const double del=1.05;
+    const double qRg = q*sqrt(Rgsquare(q,L,b));
+    const double x = pow(qRg, -1.0/miu);
+    //calculating the derivative to decide on the corection (cutoff) term?
+    // I have modified this from WRs original code
+    const double qdel = (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q);
+    const double C_star2 =(qdel >= 0.0) ? 0.0 : 1.0;
+    yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) +
+         C_star2*w_WR(qRg)*
+         x*(C1 + x*(C2 + x*C3));
+    return (yy);
+}
+double Sk_WR(double q, double L, double b)
+{
+    //
+    const double p1 = 4.12;
+    const double p2 = 4.42;
+    const double p1short = 5.36;
+    const double p2short = 5.62;
+    const double q0 = 3.1;
+    double q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
+    double Sexvmodify, ans;
+    const double C =(L/b > 10.0) ? 3.06/pow((L/b),0.44) : 1.0;
+    if( L > 4*b ) { // Longer Chains
+       if (q*b <= 3.1) {   //Modified by Yun on Oct. 15,
+         Sexvmodify = Sexvnew(q, L, b);
+         ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) -
+            (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L);
+       } else { //q(i)*b > 3.1
+         ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) +
+               a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + M_PI/(q*L);
+       }
+    } else { //L <= 4*b Shorter Chains
+       if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
+         if (q*b<=0.01) {
+            ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
+         } else {
+            ans = Sdebye1(q,L,b);
+         }
+       } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
+         ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) +
+               a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) +
+               M_PI/(q*L);
+       }
+    const double qdel = (Sexv(q*del,L,b) - Sexv_orig)/(q*(del - 1.0));
+    if (qdel < 0) {
+        //printf("branch A1-%d q=%g L=%g b=%g\n", C==1.0, q, L, b);
+        return t9 + Sexv_orig;
+    } else {
+        //printf("branch A2-%d q=%g L=%g b=%g\n", C==1.0, q, L, b);
+        const double w = w_WR(qr);
+        const double t10 = Sdebye(qr*qr)*(1.0 - w);
+        return t9 + t10;
+    }
+  return(ans);
+}
+}
+static double
+Sk_WR(double q, double L, double b)
+{
+    const double Rg_short = sqrt(Rgsquareshort(L, b));
+    double q0short = fmax(1.9/Rg_short, 3.0);
+    double ans;
+    if( L > 4*b ) { // L > 4*b : Longer Chains
+        if (q*b <= 3.1) {
+            ans = Sexv_new(q, L, b);
+        } else { //q(i)*b > 3.1
+            ans = a_long(q, L, b /*, p1, p2, q0*/);
+        }
+    } else { // L <= 4*b : Shorter Chains
+        if (q*b <= q0short) { // q*b <= fmax(1.9/Rg_short, 3)
+            //printf("branch C-%d q=%g L=%g b=%g\n", square(q*Rg_short)<DEBYE_CUTOFF, q, L, b);
+            // Note that q0short is usually 3, but it will be greater than 3
+            // small enough b, depending on the L/b ratio:
+            //     L/b == 1 => b < 2.37
+            //     L/b == 2 => b < 1.36
+            //     L/b == 3 => b < 1.00
+            //     L/b == 4 => b < 0.816
+            // 2017-10-01 pkienzle: moved low q approximation into Sdebye()
+            ans = Sdebye(square(q*Rg_short));
+        } else {  // q*b > max(1.9/Rg_short, 3)
+            //printf("branch D q=%g L=%g b=%g\n", q, L, b);
+            ans = a_short(q, L, b /*, p1short, p2short*/, q0short);
+        }
+    }
+    return ans;
+}

sasmodels/models/line.py

-                      r3330bb4
+                      rc63a7c8
 .. math::
     I(q) = \text{scale} (I(qx) \cdot I(qy)) + \text{background}
 …
 Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values
+def random():
+    import numpy as np
+    scale = 10**np.random.uniform(0, 3)
+    slope = np.random.uniform(-1, 1)*1e2
+    offset = 1e-5 + (0 if slope > 0 else -slope)
+    intercept = 10**np.random.uniform(0, 1) + offset
+    pars = dict(
+        #background=0,
+        scale=scale,
+        slope=slope,
+        intercept=intercept,
+    )
+    return pars
 tests = [

sasmodels/models/linear_pearls.py

-                      rc3ccaec
+                      r8f04da4
 source = ["lib/sas_3j1x_x.c", "linear_pearls.c"]
+demo = dict(scale=1.0, background=0.0,
+            radius=80.0,
+            edge_sep=350.0,
+            num_pearls=3,
+            sld=1.0,
+            sld_solvent=6.3)
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(1, 3) # 1 - 1000
+    edge_sep = 10**np.random.uniform(0, 3)  # 1 - 1000
+    num_pearls = np.round(10**np.random.uniform(0.3, 3)) # 2 - 1000
+    pars = dict(
+        radius=radius,
+        edge_sep=edge_sep,
+        num_pearls=num_pearls,
+    )
+    return pars
 """

sasmodels/models/lorentz.py

-                      r2c74c11
+                      r404ebbd
 description = """
 Model that evaluates a Lorentz (Ornstein-Zernicke) model.
 I(q) = scale/( 1 + (q*L)^2 ) + bkd
 The model has three parameters:
     length     =  screening Length
     scale  =  scale factor
     background    =  incoherent background
+The model has three parameters:
+    length = screening Length
+    scale = scale factor
+    background = incoherent background
 """
 category = "shape-independent"
 …
 """
+def random():
+    import numpy as np
+    pars = dict(
+        #background=0,
+        scale=10**np.random.uniform(1, 4),
+        cor_length=10**np.random.uniform(0, 3),
+    )
+    return pars
 # parameters for demo
 demo = dict(scale=1.0, background=0.0, cor_length=50.0)

sasmodels/models/mass_fractal.py

-                      r925ad6e
+                      r4553dae
 .. math::
     scale = scale\_factor \times NV^2(\rho_{particle} - \rho_{solvent})^2
+    scale = scale\_factor \times NV^2(\rho_\text{particle} - \rho_\text{solvent})^2
 .. math::
 …
 where $R$ is the radius of the building block, $D_m$ is the **mass** fractal
 dimension, | \zeta\|  is the cut-off length, $\rho_{solvent}$ is the scattering
 length density of the solvent,
 and $\rho_{particle}$ is the scattering length density of particles.
+dimension, $\zeta$  is the cut-off length, $\rho_\text{solvent}$ is the scattering
+length density of the solvent, and $\rho_\text{particle}$ is the scattering
+length density of particles.
 .. note::
     The mass fractal dimension ( $D_m$ ) is only
     valid if $0 < mass\_dim < 6$. It is also only valid over a limited
+    valid if $1 < mass\_dim < 6$. It is also only valid over a limited
     $q$ range (see the reference for details).
 …
 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "mass_fractal.c"]
+def random():
+    import numpy as np
+    radius = 10**np.random.uniform(0.7, 4)
+    cutoff_length = 10**np.random.uniform(0.7, 2)*radius
+    # TODO: fractal dimension should range from 1 to 6
+    fractal_dim_mass = 2*np.random.beta(3, 4) + 1
+    Vf = 10**np.random.uniform(-4, -1)
+    pars = dict(
+        #background=0,
+        scale=1, #1e5*Vf/radius**(fractal_dim_mass),
+        radius=radius,
+        cutoff_length=cutoff_length,
+        fractal_dim_mass=fractal_dim_mass,
+    )
+    return pars
 demo = dict(scale=1, background=0,
             radius=10.0,

sasmodels/models/mass_surface_fractal.py

-                      r6d96b66
+                      rca04add
 .. math::
+    I(q) = scale \times P(q) + background
+    I(q) = scale \times P(q) + background \\
     P(q) = \left\{ \left[ 1+(q^2a)\right]^{D_m/2} \times
                    \left[ 1+(q^2b)\right]^{(6-D_s-D_m)/2}
+           \right\}^{-1}
+    a = R_{g}^2/(3D_m/2)
+    b = r_{g}^2/[-3(D_s+D_m-6)/2]
+           \right\}^{-1} \\
+    a = R_{g}^2/(3D_m/2) \\
+    b = r_{g}^2/[-3(D_s+D_m-6)/2] \\
     scale = scale\_factor \times NV^2 (\rho_{particle} - \rho_{solvent})^2
 …
 source = ["mass_surface_fractal.c"]
+def random():
+    import numpy as np
+    fractal_dim = np.random.uniform(0, 6)
+    surface_portion = np.random.uniform(0, 1)
+    fractal_dim_surf = fractal_dim*surface_portion
+    fractal_dim_mass = fractal_dim - fractal_dim_surf
+    rg_cluster = 10**np.random.uniform(1, 5)
+    rg_primary = rg_cluster*10**np.random.uniform(-4, -1)
+    scale = 10**np.random.uniform(2, 5)
+    pars = dict(
+        #background=0,
+        scale=scale,
+        fractal_dim_mass=fractal_dim_mass,
+        fractal_dim_surf=fractal_dim_surf,
+        rg_cluster=rg_cluster,
+        rg_primary=rg_primary,
+    )
+    return pars
 demo = dict(scale=1, background=0,
             fractal_dim_mass=1.8,

sasmodels/models/mono_gauss_coil.py

-                      r3330bb4
+                      rca04add
      I_0 &= \phi_\text{poly} \cdot V
+            \cdot (\rho_\text{poly} - \rho_\text{solv})^2
+     P(q) &= 2 [\exp(-Z) + Z - 1] / Z^2
+     Z &= (q R_g)^2
+            \cdot (\rho_\text{poly} - \rho_\text{solv})^2 \\
+     P(q) &= 2 [\exp(-Z) + Z - 1] / Z^2 \\
+     Z &= (q R_g)^2 \\
      V &= M / (N_A \delta)
 …
 description = """
     Evaluates the scattering from
+    Evaluates the scattering from
     monodisperse polymer chains.
     """
 …
 Iq.vectorized = True # Iq accepts an array of q values
+def random():
+    import numpy as np
+    rg = 10**np.random.uniform(0, 4)
+    #rg = 1e3
+    pars = dict(
+        #scale=1, background=0,
+        i_zero=1e7, # i_zero is a simple scale
+        rg=rg,
+    )
+    return pars
 demo = dict(scale=1.0, i_zero=70.0, rg=75.0, background=0.0)

sasmodels/models/multilayer_vesicle.py

-                      r870a2f4
+                      r64eecf7
 .. math::
      r_i &= r_c + (i-1)(t_s + t_w) && \text{ solvent radius before shell } i \\
      R_i &= r_i + t_s && \text{ shell radius for shell } i
+     r_i &= r_c + (i-1)(t_s + t_w) \text{ solvent radius before shell } i \\
+     R_i &= r_i + t_s \text{ shell radius for shell } i
 $\phi$ is the volume fraction of particles, $V(r)$ is the volume of a sphere
 …
     return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent
+demo = dict(scale=1, background=0,
+            volfraction=0.05,
+            radius=60.0,
+            thick_shell=10.0,
+            thick_solvent=10.0,
+            sld_solvent=6.4,
+            sld=0.4,
+            n_shells=2.0)
+def random():
+    import numpy as np
+    volfraction = 10**np.random.uniform(-3, -0.5)  # scale from 0.1% to 30%
+    radius = 10**np.random.uniform(0, 2.5) # core less than 300 A
+    total_thick = 10**np.random.uniform(2, 4) # up to 10000 A of shells
+    # random number of shells, with shell+solvent thickness > 10 A
+    n_shells = int(10**np.random.uniform(0, np.log10(total_thick)-1)+0.5)
+    # split total shell thickness with preference for shell over solvent;
+    # make sure that shell thickness is at least 1 A
+    one_thick = total_thick/n_shells
+    thick_solvent = 10**np.random.uniform(-2, 0)*(one_thick - 1)
+    thick_shell = one_thick - thick_solvent
+    pars = dict(
+        scale=1,
+        volfraction=volfraction,
+        radius=radius,�

Changeset 41917b8 in sasmodels

Legend: