Changeset bd91e8f – SasView

.travis.yml

-                      r5c36bf1
+                      rf3767c2
   - os: linux
     env:
     - PY=3.6
+    - PY=3.7
   - os: osx
     language: generic
 …
     language: generic
     env:
     - PY=3.5
+    - PY=3.7
 branches:
   only:

doc/guide/index.rst

-                      rda5536f
+                      rbc69321
    pd/polydispersity.rst
    resolution.rst
+   plugin.rst
+   fitting_sq.rst
    magnetism/magnetism.rst
    orientation/orientation.rst
    sesans/sans_to_sesans.rst
    sesans/sesans_fitting.rst
-   plugin.rst
    scripting.rst
    refs.rst

doc/guide/plugin.rst

-                      r9150036
+                      re15a822
 **Note: The order of the parameters in the definition will be the order of the
 parameters in the user interface and the order of the parameters in Fq(), Iq(),
 Iqac(), Iqabc(), form_volume() and shell_volume().
+Iqac(), Iqabc(), radius_effective(), form_volume() and shell_volume().
 And** *scale* **and** *background* **parameters are implicit to all models,
 so they do not need to be included in the parameter table.**
 …
 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::
+you should round to the nearest integer.  In C code, you can do this using:
+.. code-block:: c
     static double
 …
 .............
+.. note::
+   Pure python models do not yet support direct computation of $<F(Q)^2>$ or
+   $<F(Q)>^2$. Neither do they support orientational distributions or magnetism
+   (use C models if these are required).
 For pure python models, define the *Iq* function::
 …
 Models should define *form_volume(par1, par2, ...)* where the parameter
 list includes the *volume* parameters in order.  This is used for a weighted
+volume normalization so that scattering is on an absolute scale.  If
+*form_volume* is not defined, then the default *form_volume = 1.0* will be
+used.
+volume normalization so that scattering is on an absolute scale.  For
+solid shapes, the *I(q)* function should use *form_volume* squared
+as its scale factor.  If *form_volume* is not defined, then the default
+*form_volume = 1.0* will be used.
 Hollow shapes, where the volume fraction of particle corresponds to the
 material in the shell rather than the volume enclosed by the shape, must
 also define a *shell_volume(par1, par2, ...)* function.  The parameters
+are the same as for *form_volume*.  The *I(q)* calculation should use
+*shell_volume* squared as its scale factor for the volume normalization.
+The structure factor calculation needs *form_volume* in order to properly
+scale the volume fraction parameter, so both functions are required for
+hollow shapes.
+Note: Pure python models do not yet support direct computation of the
+average of $F(q)$ and $F^2(q)$.
+are the same as for *form_volume*.  Here the *I(q)* function should use
+*shell_volume* squared instead of *form_volume* squared so that the scale
+parameter corresponds to the volume fraction of material in the sample.
+The structure factor calculation needs the volume fraction of the filled
+shapes for its calculation, so the volume fraction parameter in the model
+is automatically scaled by *form_volume/shell_volume* prior to calling the
+structure factor.
 Embedded C Models
 …
     """
+This expands into the equivalent C code::
+This expands into the equivalent C code:
+.. code-block:: c
     double Iq(double q, double par1, double par2, ...);
 …
 includes only the volume parameters.
 *form_volume* defines the volume of the shell for hollow shapes. As in
+*shell_volume* defines the volume of the shell for hollow shapes. As in
 python models, it includes only the volume parameters.
 …
 Rather than returning NAN from Iq, you must define the *INVALID(v)*.  The
 *v* parameter lets you access all the parameters in the model using
+*v.par1*, *v.par2*, etc. For example::
+*v.par1*, *v.par2*, etc. For example:
+.. code-block:: c
     #define INVALID(v) (v.bell_radius < v.radius)
 …
 Instead of defining the *Iq* function, models can define *Fq* as
+something like::
+something like:
+.. code-block:: c
     double Fq(double q, double *F1, double *F2, double par1, double par2, ...);
 …
 laboratory frame and beam travelling along $-z$.
+The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where
+The oriented C model (oriented pure Python models are not supported)
+is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where
 *par1*, etc. are the parameters to the model.  If the shape is rotationally
 symmetric about *c* then *psi* is not needed, and the model is called
 …
 Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the
+numerical integration is then::
+numerical integration is then:
+.. code-block:: c
     double outer_sum = 0.0;
 …
 to compute the proper magnetism and orientation, which you can implement
 using *Iqxy(qx, qy, par1, par2, ...)*.
+**Note: Magnetism is not supported in pure Python models.**
 Special Functions
 …
   memory, and wrong answers computed. The conclusion from a very long and
   strange debugging session was that any arrays that you declare in your
+  model should be a multiple of four. For example::
+  model should be a multiple of four. For example:
+  .. code-block:: c
       double Iq(q, p1, p2, ...)
 …
 structure factor is the *hardsphere* interaction, which
 uses the effective radius of the form factor as an input to the structure
+factor model.  The effective radius is the average radius of the
+form averaged over all the polydispersity values.
+::
+    def ER(radius, thickness):
+        """Effective radius of a core-shell sphere."""
+        return radius + thickness
+Now consider the *core_shell_sphere*, which has a simple effective radius
+equal to the radius of the core plus the thickness of the shell, as
+shown above. Given polydispersity over *(r1, r2, ..., rm)* in radius and
+*(t1, t2, ..., tn)* in thickness, *ER* is called with a mesh
+grid covering all possible combinations of radius and thickness.
+That is, *radius* is *(r1, r2, ..., rm, r1, r2, ..., rm, ...)*
+and *thickness* is *(t1, t1, ... t1, t2, t2, ..., t2, ...)*.
+The *ER* function returns one effective radius for each combination.
+The effective radius calculator weights each of these according to
+the polydispersity distributions and calls the structure factor
+with the average *ER*.
+::
+    def VR(radius, thickness):
+        """Sphere and shell volumes for a core-shell sphere."""
+        whole = 4.0/3.0 * pi * (radius + thickness)**3
+        core = 4.0/3.0 * pi * radius**3
+        return whole, whole - core
+Core-shell type models have an additional volume ratio which scales
+the structure factor.  The *VR* function returns the volume of
+the whole sphere and the volume of the shell. Like *ER*, there is
+one return value for each point in the mesh grid.
+*NOTE: we may be removing or modifying this feature soon. As of the
+time of writing, core-shell sphere returns (1., 1.) for VR, giving a volume
+ratio of 1.0.*
+factor model.  The effective radius is the weighted average over all
+values of the shape in polydisperse systems.
+There can be many notions of effective radius, depending on the shape.  For
+a sphere it is clearly just the radius, but for an ellipsoid of revolution
+we provide average curvature, equivalent sphere radius, minimum radius and
+maximum radius.  These options are listed as *radius_effective_modes* in
+the python model defintion, and must be computed by the *radius_effective*
+function which takes the *radius_effective_mode* parameter as an integer,
+along with the various model parameters.  Unlike normal C/Python arrays,
+the first mode is 1, the second is 2, etc.  Mode 0 indicates that the
+effective radius from the shape is to be ignored in favour of the the
+effective radius parameter in the structure factor model.
+Consider the core-shell sphere, which defines the following effective radius
+modes in the python model::
+    radius_effective_modes = [
+        "outer radius",
+        "core radius",
+    ]
+and the following function in the C-file for the model:
+.. code-block:: c
+    static double
+    radius_effective(int mode, double radius, double thickness)
+    {
+        switch (mode) {
+            case 0: return radius + thickness;
+            case 1: return radius;
+            default: return 0.;
+        }
+    }
+    static double
+    form_volume(double radius, double thickness)
+    {
+        return M_4PI_3 * cube(radius + thickness);
+    }
+Given polydispersity over *(r1, r2, ..., rm)* in radius and *(t1, t2, ..., tn)*
+in thickness, *radius_effective* is called over a mesh grid covering all
+possible combinations of radius and thickness, with a single *(ri, tj)* pair
+in each call. The weights each of these results according to the
+polydispersity distributions and calls the structure factor with the average
+effective radius.  Similarly, for *form_volume*.
+Hollow models have an additional volume ratio which is needed to scale the
+structure factor.  The structure factor uses the volume fraction of the filled
+particles as part of its density estimate, but the scale factor for the
+scattering intensity (as non-solvent volume fraction / volume) is determined
+by the shell volume only.  Therefore the *shell_volume* function is
+needed to compute the form:shell volume ratio, which then scales the
+*volfraction* parameter prior to calling the structure factor calculator.
+In the case of a hollow sphere, this would be:
+.. code-block:: c
+    static double
+    shell_volume(double radius, double thickness)
+    {
+        double whole = M_4PI_3 * cube(radius + thickness);
+        double core = M_4PI_3 * cube(radius);
+        return whole - core;
+    }
+If *shell_volume* is not present, then *form_volume* and *shell_volume* are
+assumed to be equal, and the shape is considered solid.
 Unit Tests
 …
 and a check that the model runs.
-If you are not using sasmodels from SasView, skip this step.
 Recommended Testing
 ...................
+**NB: For now, this more detailed testing is only possible if you have a
+SasView build environment available!**
 If the model compiles and runs, you can next run the unit tests that
 …
 | 2016-10-25 Steve King
 | 2017-05-07 Paul Kienzle - Moved from sasview to sasmodels docs
+| 2019-03-28 Paul Kienzle - Update docs for radius_effective and shell_volume

doc/guide/resolution.rst

-                      r0db85af
+                      rdb1d9d5
 .. sm_help.rst
 .. This is a port of the original SasView html help file to ReSTructured text
 .. by S King, ISIS, during SasView CodeCamp-III in Feb 2015.
+.. resolution.rst
+.. This is a port of the original SasView html help file sm_help to ReSTructured
+.. text by S King, ISIS, during SasView CodeCamp-III in Feb 2015.
 …
 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*. The Sasmodels component of SasView
+does the latter.
 Both smearing and desmearing rely on functions to describe the resolution
 …
 for the instrument and stored with the data file.  If not, they will need to
 be set manually before fitting.
+.. note::
+    Problems may be encountered if the data set loaded by SasView is a
+    concatenation of SANS data from several detector distances where, of
+    course, the worst Q resolution is next to the beam stop at each detector
+    distance. (This will also be noticeable in the residuals plot where
+    there will be poor overlap). SasView sensibly orders all the input
+    data points by increasing Q for nicer-looking plots, however, the dQ
+    data can then vary considerably from point to point. If 'Use dQ data'
+    smearing is selected then spikes may appear in the model fits, whereas
+    if 'None' or 'Custom Pinhole Smear' are selected the fits look normal.
+    In such instances, possible solutions are to simply remove the data
+    with poor Q resolution from the shorter detector distances, or to fit
+    the data from different detector distances simultaneously.

explore/asymint.py

-                      rc462169
+                      r66dbbfb
 import os, sys
 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
+import warnings
 import numpy as np
 …
 import sasmodels.special as sp
-# Need to parse shape early since it determines the kernel function
-# that will be used for the various integrators
-shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1]
-Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2]
 DTYPE = 'd'
 …
     return norm, Fq
+if shape == 'sphere':
+    RADIUS = 50  # integer for the sake of mpf
+    NORM, KERNEL = make_sphere(radius=RADIUS)
+    NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv)
+elif shape == 'cylinder':
+    #RADIUS, LENGTH = 10, 100000
+    RADIUS, LENGTH = 10, 300  # integer for the sake of mpf
+    NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH)
+    NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv)
+elif shape == 'triaxial_ellipsoid':
+    #A, B, C = 4450, 14000, 47
+    A, B, C = 445, 140, 47  # integer for the sake of mpf
+    NORM, KERNEL = make_triaxial_ellipsoid(A, B, C)
+    NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv)
+elif shape == 'parallelepiped':
+    #A, B, C = 4450, 14000, 47
+    A, B, C = 445, 140, 47  # integer for the sake of mpf
+    NORM, KERNEL = make_parallelepiped(A, B, C)
+    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
+elif shape == 'core_shell_parallelepiped':
+    #A, B, C = 4450, 14000, 47
+    #A, B, C = 445, 140, 47  # integer for the sake of mpf
+    A, B, C = 114, 1380, 6800
+    DA, DB, DC = 21, 58, 2300
+    SLDA, SLDB, SLDC = "5", "-0.3", "11.5"
+    ## default parameters from sasmodels
+    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2
+    ## swap A-B-C to C-B-A
+    #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA
+    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3
+    #SLD_SOLVENT,CONTRAST = 0, 4
+    if 1: # C shortest
+        B, C = C, B
+        DB, DC = DC, DB
+        SLDB, SLDC = SLDC, SLDB
+    elif 0: # C longest
+        A, C = C, A
+        DA, DC = DC, DA
+        SLDA, SLDC = SLDC, SLDA
+    #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
+    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
+    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
+elif shape == 'paracrystal':
+    LATTICE = 'bcc'
+    #LATTICE = 'fcc'
+    #LATTICE = 'sc'
+    DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string
+    RADIUS = 40  # integer for the sake of mpf
+    NORM, KERNEL = make_paracrystal(
+        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE)
+    NORM_MP, KERNEL_MP = make_paracrystal(
+        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv)
+else:
+    raise ValueError("Unknown shape %r"%shape)
+NORM = 1.0  # type: float
+KERNEL = None  # type: CALLABLE[[ndarray, ndarray, ndarray], ndarray]
+NORM_MP = 1  # type: mpf
+KERNEL = None  # type: CALLABLE[[mpf, mpf, mpf], mpf]
+SHAPES = [
+    'sphere',
+    'cylinder',
+    'triaxial_ellipsoid',
+    'parallelepiped',
+    'core_shell_parallelepiped',
+    'fcc_paracrystal',
+    'bcc_paracrystal',
+    'sc_paracrystal',
+]
+def build_shape(shape, **pars):
+    global NORM, KERNEL
+    global NORM_MP, KERNEL_MP
+    # Note: using integer or string defaults for the sake of mpf
+    if shape == 'sphere':
+        RADIUS = pars.get('radius', 50)
+        NORM, KERNEL = make_sphere(radius=RADIUS)
+        NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv)
+    elif shape == 'cylinder':
+        #RADIUS, LENGTH = 10, 100000
+        RADIUS = pars.get('radius', 10)
+        LENGTH = pars.get('radius', 300)
+        NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH)
+        NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv)
+    elif shape == 'triaxial_ellipsoid':
+        #A, B, C = 4450, 14000, 47
+        A = pars.get('a', 445)
+        B = pars.get('b', 140)
+        C = pars.get('c', 47)
+        NORM, KERNEL = make_triaxial_ellipsoid(A, B, C)
+        NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv)
+    elif shape == 'parallelepiped':
+        #A, B, C = 4450, 14000, 47
+        A = pars.get('a', 445)
+        B = pars.get('b', 140)
+        C = pars.get('c', 47)
+        NORM, KERNEL = make_parallelepiped(A, B, C)
+        NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
+    elif shape == 'core_shell_parallelepiped':
+        #A, B, C = 4450, 14000, 47
+        #A, B, C = 445, 140, 47  # integer for the sake of mpf
+        A = pars.get('a', 114)
+        B = pars.get('b', 1380)
+        C = pars.get('c', 6800)
+        DA = pars.get('da', 21)
+        DB = pars.get('db', 58)
+        DC = pars.get('dc', 2300)
+        SLDA = pars.get('slda', "5")
+        SLDB = pars.get('sldb', "-0.3")
+        SLDC = pars.get('sldc', "11.5")
+        ## default parameters from sasmodels
+        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2
+        ## swap A-B-C to C-B-A
+        #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA
+        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3
+        #SLD_SOLVENT,CONTRAST = 0, 4
+        if 1: # C shortest
+            B, C = C, B
+            DB, DC = DC, DB
+            SLDB, SLDC = SLDC, SLDB
+        elif 0: # C longest
+            A, C = C, A
+            DA, DC = DC, DA
+            SLDA, SLDC = SLDC, SLDA
+        #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
+        NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
+        NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
+    elif shape.endswith('paracrystal'):
+        LATTICE, _ = shape.split('_')
+        DNN = pars.get('dnn', 220)
+        D_FACTOR = pars.get('d_factor', '0.06')
+        RADIUS = pars.get('radius', 40)
+        NORM, KERNEL = make_paracrystal(
+            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE)
+        NORM_MP, KERNEL_MP = make_paracrystal(
+            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv)
+    else:
+        raise ValueError("Unknown shape %r"%shape)
 # Note: hardcoded in mp_quad
 …
 # 2D integration functions
 def mp_quad_2d(q, shape):
+def mp_quad_2d(q):
     evals = [0]
     def integrand(theta, phi):
 …
     print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
     print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
+def quadpy_method(q, rule):
+    """
+    Use *rule*="name:index" where name and index are chosen from below.
+    Available rule names and the corresponding indices::
+        AlbrechtCollatz: [1-5]
+        BazantOh: 9, 11, 13
+        HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3],
+, 31, 33, 35, 37, 39-[1-2]
+        FliegeMaier: 4, 9, 16, 25
+        Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35,
+, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131
+        McLaren: [1-10]
+        Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3],
+            U3 11-[1-3], U3 14-1
+    """
+    try:
+        import quadpy
+    except ImportError:
+        warnings.warn("use 'pip install quadpy' to enable quadpy.sphere tests")
+        return
+    from quadpy.sphere import (AlbrechtCollatz, BazantOh, HeoXu,
+        FliegeMaier, Lebedev, McLaren, Stroud, integrate_spherical)
+    RULES = {
+        'AlbrechtCollatz': AlbrechtCollatz,
+        'BazantOh': BazantOh,
+        'HeoXu': HeoXu,
+        'FliegeMaier': FliegeMaier,
+        'Lebedev': Lebedev,
+        'McLaren': McLaren,
+        'Stroud': Stroud,
+    }
+    int_index = 'AlbrechtCollatz', 'McLaren'
+    rule_name, rule_index = rule.split(':')
+    index = int(rule_index) if rule_name in int_index else rule_index
+    rule_obj = RULES[rule_name](index)
+    fn = lambda azimuthal, polar: kernel_2d(q=q, theta=polar, phi=azimuthal)
+    Iq = integrate_spherical(fn, rule=rule_obj)/(4*pi)
+    print("%s degree=%d points=%s => %.15g"
+          % (rule, rule_obj.degree, len(rule_obj.points), Iq))
 def plot_2d(q, n=300):
 …
     pylab.show()
+def main(Qstr):
+    Q = float(Qstr)
+    if shape == 'sphere':
+def main():
+    import argparse
+    parser = argparse.ArgumentParser(
+        description="asymmetric integration explorer",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+        )
+    parser.add_argument('-s', '--shape', choices=SHAPES,
+                        default='parallelepiped',
+                        help='oriented shape')
+    parser.add_argument('-q', '--q_value', type=str, default='0.005',
+                        help='Q value to evaluate')
+    parser.add_argument('pars', type=str, nargs='*', default=[],
+                        help='p=val for p in shape parameters')
+    opts = parser.parse_args()
+    pars = {k: v for par in opts.pars for k, v in [par.split('=')]}
+    build_shape(opts.shape, **pars)
+    Q = float(opts.q_value)
+    if opts.shape == 'sphere':
         print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2)
+    print("gauss-20", *gauss_quad_2d(Q, n=20))
+    print("gauss-76", *gauss_quad_2d(Q, n=76))
+    print("gauss-150", *gauss_quad_2d(Q, n=150))
+    print("gauss-500", *gauss_quad_2d(Q, n=500))
+    print("gauss-1025", *gauss_quad_2d(Q, n=1025))
+    print("gauss-2049", *gauss_quad_2d(Q, n=2049))
+    print("gauss-20 usub", *gauss_quad_usub(Q, n=20))
+    print("gauss-76 usub", *gauss_quad_usub(Q, n=76))
+    print("gauss-150 usub", *gauss_quad_usub(Q, n=150))
+    # Methods from quadpy, if quadpy is available
+    #  AlbrechtCollatz: [1-5]
+    #  BazantOh: 9, 11, 13
+    #  HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3],
+    #     29, 31, 33, 35, 37, 39-[1-2]
+    #  FliegeMaier: 4, 9, 16, 25
+    #  Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35,
+    #     41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131
+    #  McLaren: [1-10]
+    #  Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3],
+    #     U3 11-[1-3], U3 14-1
+    quadpy_method(Q, "AlbrechtCollatz:5")
+    quadpy_method(Q, "HeoXu:39-2")
+    quadpy_method(Q, "FliegeMaier:25")
+    quadpy_method(Q, "Lebedev:19")
+    quadpy_method(Q, "Lebedev:131")
+    quadpy_method(Q, "McLaren:10")
+    quadpy_method(Q, "Stroud:U3 14-1")
+    print("gauss-20 points=%d => %.15g" % gauss_quad_2d(Q, n=20))
+    print("gauss-76 points=%d => %.15g" % gauss_quad_2d(Q, n=76))
+    print("gauss-150 points=%d => %.15g" % gauss_quad_2d(Q, n=150))
+    print("gauss-500 points=%d => %.15g" % gauss_quad_2d(Q, n=500))
+    print("gauss-1025 points=%d => %.15g" % gauss_quad_2d(Q, n=1025))
+    print("gauss-2049 points=%d => %.15g" % gauss_quad_2d(Q, n=2049))
+    print("gauss-20 usub points=%d => %.15g" % gauss_quad_usub(Q, n=20))
+    print("gauss-76 usub points=%d => %.15g" % gauss_quad_usub(Q, n=76))
+    print("gauss-150 usub points=%d => %.15g" % gauss_quad_usub(Q, n=150))
     #gridded_2d(Q, n=2**8+1)
     gridded_2d(Q, n=2**10+1)
     #gridded_2d(Q, n=2**12+1)
     #gridded_2d(Q, n=2**15+1)
+    if shape not in ('paracrystal', 'core_shell_parallelepiped'):
+        # adaptive forms on models for which the calculations are fast enough
+    # adaptive forms on models for which the calculations are fast enough
+    SLOW_SHAPES = {
+        'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal',
+        'core_shell_parallelepiped',
+    }
+    if opts.shape not in SLOW_SHAPES:
         print("dblquad", *scipy_dblquad_2d(Q))
         print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
         print("romberg", *scipy_romberg_2d(Q))
         with mp.workprec(100):
             print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape))
+            print("mpmath", *mp_quad_2d(mp.mpf(opts.q_value)))
     plot_2d(Q, n=200)
 if __name__ == "__main__":
     main(Qstr)
+    main()

explore/beta/sasfit_compare.py

r119073a	ra34b811
408	408	#radius_effective=12.59921049894873,
409	409	)
410		target = sasmodels_theory(q, model, ~~effective_radius~~_mode=0, structure_factor_mode=1, **pars)
	410	target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars)
411	411	actual = ellipsoid_pe(q, norm='sasview', **pars)
412	412	title = " ".join(("sasmodels", model, pd_type))

explore/beta/sasfit_compare_new.py

-                      r2586ab72
+                      ra34b811
 #Used to calculate F(q) for the cylinder, sphere, ellipsoid models
 # RKH adding vesicle and hollow_cylinder to test sasview special cases of ER and VR
 # BEWARE there are issues here if you run this from python3 (i.e. qt5 version of sasview), so do "source activate sasview"
+# There were issues here from python3 (i.e. qt5 version of sasview),fixed by Paul K (else do "source activate sasview")
 def sas_sinx_x(x):
     with np.errstate(all='ignore'):
 …
     if radius_effective is None:
         radius_effective = radius
-        print("radius_effective  ",radius_effective)
         average_volume = total_volume/total_weight
     if norm == 'sasfit':
 …
         F2 *= 1e-12
         IQD = F2/average_volume*1e8*volfraction
+    print("in sphere_r : radius_effective  ",radius_effective,"  volfraction",volfraction)
     SQ = hardsphere_simple(q, radius_effective, volfraction)
     beta = F1**2/F2
 …
     if radius_effective is None:
         radius_effective = radius_eff/total_weight
         print("radius_effective  ",radius_effective)
+        print("in vesicle_pe : radius_effective  ",radius_effective)
     if norm == 'sasfit':
         IQD = F2
 …
         F2 *= 1e-12
         IQD = F2/average_volume*1e8*volfraction
     # RKH SORT THIS OUT WHEN HAVE NEW MODEL
     vfff=volfraction*(30./20.)**3
     print("radius_effective, fudged volfaction for S(Q) ",radius_effective,vfff)
+    # RKH SORT THIS OUT WHEN HAVE NEW MODEL  Vtot/Vshell = (R+T)^3/ ( (R+T)^3 -R^3 )
+    vfff=volfraction*( 1.0/(1.0 - ( (radius/(radius+thickness))**3 )))
+    print("in vesicle_pe : radius_effective, fudged volfraction for S(Q) ",radius_effective,vfff)
     SQ = hardsphere_simple(q, radius_effective, vfff)
     beta = F1**2/F2
 …
+#
 #polydispersity for hollow_cylinder
 def hollow_cylinder_pe(q,radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.15,
         radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", thickness_pd_type="gaussian",
         radius_effective=None, background=0, scale=1,
+def hollow_cylinder_pe(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1, volfraction=0.15,
+        radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", length_pd_type="gaussian",
+        thickness_pd_type="gaussian", radius_effective=None, background=0, scale=1,
                  norm='sasview'):
     if norm not in ['sasview', 'sasfit', 'yun']:
 …
         T_val, T_prob = schulz_distribution(thickness, thickness_pd*thickness, 0, inf)
     if length_pd_type == 'gaussian':
         T_val, T_prob = gaussian_distribution(length, length_pd*length, 0, inf)
+        L_val, L_prob = gaussian_distribution(length, length_pd*length, 0, inf)
     elif length_pd_type == 'schulz':
         L_val, L_prob = schulz_distribution(length, length_pd*length, 0, inf)
 …
     if radius_effective is None:
         radius_effective = radius_eff/total_weight
         print("radius_effective  ",radius_effective)
+        print("in hollow_cylinder : radius_effective  ",radius_effective)
     if norm == 'sasfit':
         IQD = F2
 …
 # RKH SORT THIS OUT WHEN HAVE NEW MODEL
     vfff=volfraction
     print("radius_effective, volfaction for S(Q) ",radius_effective,vfff)
+    print("in hollow_cylinder : radius_effective, volfaction for S(Q) ",radius_effective,vfff)
     SQ = hardsphere_simple(q, radius_effective, vfff)
     beta = F1**2/F2
 …
     #Iq = None#= Sq = None
     r = dict(I._kernel.results())
     #print(r)
+    print(r)
     return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq)
 …
         radius_pd=.1, radius_pd_type=pd_type,
         volfraction=0.15,
+        radius_effective=12.59921049894873,  # equivalent average sphere radius
+        radius_effective=20  # equivalent average sphere radius
+#       NOTE sasview computes its own radius_effective in "target" (the print(r) at end of sasmodels_theory will reveal its value),
+#       this one is only used locally for "actual"
+        )
     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
+    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars)
     actual = sphere_r(q, norm='sasview', **pars)
     title = " ".join(("sasmodels", model, pd_type))
 …
     q = np.logspace(-5, 0, 250)
     model = 'vesicle'
     print("F2 seems OK, but big issues with S(Q), wo work in progress")
 # NOTE parameers for vesicle model are soon to be changed to remove "volfraction"
+    print("F2 seems OK, but big issues with S(Q), so work in progress")
+# NOTE parameters for vesicle model are soon to be changed to remove "volfraction" (though still there in 5.0)
     pars = dict(
         radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03,
         background=0,
+        radius_pd=0.1, thickness_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type
+        #,radius_effective=12.59921049894873,  # equivalent average sphere radius
+        )
+    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
+        radius_pd=0.01, thickness_pd=0.01, radius_pd_type=pd_type, thickness_pd_type=pd_type,
+        radius_effective=30. )
+        # equivalent average sphere radius for local "actual" to match what sasview uses, use None to compute average outer radius here,
+    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars)
     actual = vesicle_pe(q, norm='sasview', **pars)
     title = " ".join(("sasmodels", model, pd_type))
 …
         radius=20, thickness=10, length=80, sld=4, sld_solvent=1,
         background=0,
         radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type
         #,radius_effective=12.59921049894873,  # equivalent average sphere radius
+        )
     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
+        radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type,
+        radius_effective=40.687)
+        # equivalent average sphere radius for local "actual" to match what sasview uses
+    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars)
     actual = hollow_cylinder_pe(q, norm='sasview', **pars)
 # RKH monodisp was OK,    actual = hollow_cylinder_theta(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1 )
 …
         volfraction=0.15,
         radius_effective=270.7543927018,
+        #radius_effective=12.59921049894873,
+# if change radius_effective to some other value, the S(Q) from sasview does not agree
+        )
     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
+    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars)
     actual = ellipsoid_pe(q, norm='sasview', **pars)
 # RKH test       actual = ellipsoid_theta(q, radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, volfraction=0.15, radius_effective=270.)
 …
+    }
     Q, IQ = load_sasfit(data_file('richard_test.txt'))
     Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
     Q, IQBD = load_sasfit(data_file('richard_test3.txt'))
+    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt'))
+    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt'))
+    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt'))
     target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
     actual = sphere_r(Q, norm="sasfit", **pars)
 …
+    }
     Q, IQ = load_sasfit(data_file('richard_test4.txt'))
     Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
     Q, IQBD = load_sasfit(data_file('richard_test6.txt'))
+    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQD.txt'))
+    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQSD.txt'))
+    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQBD.txt'))
     target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
     actual = ellipsoid_pe(Q, norm="sasfit", **pars)

extra/build_linux.sh

rdb8756e	r0d362b7
14	14	# Run tests
15	15	STATUS=0
16		python -m sasmodels.model_test ~~opencl_and_dll~~ all \|\| STATUS=$?
	16	python -m sasmodels.model_test -v all \|\| STATUS=$?
17	17	python -m sasmodels.resolution \|\| STATUS=$?
18	18

extra/pylint.rc

-                      r6e68289
+                      rdf20715
     locally-enabled,
     locally-disabled,
+    #no-else-return,
+    #invalid-name,
+    fixme,
+    too-many-lines,
+    too-many-locals,
+    too-many-branches,
+    too-many-statements,
+    #too-many-instance-attributes,
+    #too-many-return-statements,
+    #no-self-use,
 [REPORTS]
 …
 good-names=_,
     input, id,
     h, i, j, k, n, p, v, w, x, y, z,
     q, qx, qy, qz, Q, Qx, Qy, Qz,
     dt, dx, dy, dz, dq,
+    h, i, j, k, n, p, u, v, w, x, y, z, r, dr,
+    q, qx, qy, qz, Q, Qx, Qy, Qz, nx, ny, nz, nq, nr,
+    dt, dx, dy, dz, dq, S, P, a, b, c,
     Iq, dIq, Iqxy, Iq_calc,
     ER, call_ER, VR, call_VR,
 …
 # Regular expression matching correct variable names
 variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{2,30}$
+variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{1,30}$
 # Naming hint for variable names
 …
 # Regular expression matching correct argument names
 argument-rgx=[a-z_][a-z0-9_]{2,30}$
+argument-rgx=[a-z_][a-z0-9_]{1,30}$
 # Naming hint for argument names
 …
 # (useful for modules/projects where namespaces are manipulated during runtime
 # and thus existing member attributes cannot be deduced by static analysis
 ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special
+ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special,pyopencl,pycuda,matplotlib.cm,ipyvolume,ipywidget
 # List of classes names for which member attributes should not be checked

sasmodels/bumps_model.py

-                      r2c4a190
+                      rb297ba9
     from .kernel import KernelModel
     from .modelinfo import ModelInfo
+    from .resolution import Resolution
     Data = Union[Data1D, Data2D]
 except ImportError:
 …
     @property
     def resolution(self):
+        # type: () -> Union[None, Resolution]
+        """
+        :class:`sasmodels.Resolution` applied to the data, if any.
+        """
         return self._resolution
     @resolution.setter
     def resolution(self, value):
+        # type: (Resolution) -> None
+        """
+        :class:`sasmodels.Resolution` applied to the data, if any.
+        """
         self._resolution = value
 …
         Save the model parameters and data into a file.
         Not Implemented.
+        Not Implemented except for sesans fits.
         """
         if self.data_type == "sesans":

sasmodels/compare.py

-                      rc1799d3
+                      rb297ba9
         set_integration_size(model_info, ngauss)
     if (dtype != "default" and not dtype.endswith('!')
+    if (dtype != "default" and not dtype.endswith('!')
             and not (kernelcl.use_opencl() or kernelcuda.use_cuda())):
         raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR)
 …
         opts['pars'] = parse_pars(opts, maxdim=maxdim)
         if opts['pars'] is None:
             return
+            return limits
         result = run_models(opts, verbose=True)
         if opts['plot']:
 …
     args = dict((k, v) for k, v in args.items()
                 if "_pd" not in k
                    and ":" not in k
                    and k not in ("background", "scale", "theta", "phi", "psi"))
+                and ":" not in k
+                and k not in ("background", "scale", "theta", "phi", "psi"))
     args = args.copy()
 …
     # work with trimmed data, not the full set
     sorted_err = np.sort(abs(err.compressed()))
     if len(sorted_err) == 0:
+    if sorted_err.size == 0:
         print(label + "  no valid values")
         return

sasmodels/convert.py

-                      rb3f4831
+                      rbd91e8f
                     translation[newid+str(k)] = oldid+str(k)
     # Remove control parameter from the result
+    if model_info.control:
+        translation[model_info.control] = "CONTROL"
+    control_pars = [p.id for p in model_info.parameters.kernel_parameters
+                    if p.is_control]
+    if control_pars:
+        control_id = control_pars[0]
+        translation[control_id] = "CONTROL"
     return translation
 …
 def revert_name(model_info):
+    """Translate model name back to the name used in SasView 3.x"""
     oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}])
     return oldname
 …
 def test_backward_forward():
+    """
+    Test conversion of model parameters from 4.x to 3.x and back.
+    """
     from .core import list_models
     L = lambda name: _check_one(name, seed=1)

sasmodels/core.py

-                      r07646b6
+                      r9562dd2
         * all: all models
         * py: python models only
+        * c: compiled models only
+        * single: models which support single precision
+        * double: models which require double precision
+        * opencl: controls if OpenCL is supperessed
+        * 1d: models which are 1D only, or 2D using abs(q)
+        * 2d: models which can be 2D
+        * magnetic: models with an sld
+        * nommagnetic: models without an sld
+        * c: c models only
+        * single: c models which support single precision
+        * double: c models which require double precision
+        * opencl: c models which run in opencl
+        * dll: c models which do not run in opencl
+        * 1d: models without orientation
+        * 2d: models with orientation
+        * magnetic: models supporting magnetic sld
+        * nommagnetic: models without magnetic parameter
     For multiple conditions, combine with plus.  For example, *c+single+2d*
 …
     info = load_model_info(name)
     pars = info.parameters.kernel_parameters
+    if kind == "py" and callable(info.Iq):
+        return True
+    elif kind == "c" and not callable(info.Iq):
+        return True
+    elif kind == "double" and not info.single:
+        return True
+    elif kind == "single" and info.single:
+        return True
+    elif kind == "opencl" and info.opencl:
+        return True
+    elif kind == "2d" and any(p.type == 'orientation' for p in pars):
+        return True
+    elif kind == "1d" and all(p.type != 'orientation' for p in pars):
+        return True
+    elif kind == "magnetic" and any(p.type == 'sld' for p in pars):
+        return True
+    elif kind == "nonmagnetic" and any(p.type != 'sld' for p in pars):
+        return True
+    # TODO: may be adding Fq to the list at some point
+    is_pure_py = callable(info.Iq)
+    if kind == "py":
+        return is_pure_py
+    elif kind == "c":
+        return not is_pure_py
+    elif kind == "double":
+        return not info.single and not is_pure_py
+    elif kind == "single":
+        return info.single and not is_pure_py
+    elif kind == "opencl":
+        return info.opencl
+    elif kind == "dll":
+        return not info.opencl and not is_pure_py
+    elif kind == "2d":
+        return any(p.type == 'orientation' for p in pars)
+    elif kind == "1d":
+        return all(p.type != 'orientation' for p in pars)
+    elif kind == "magnetic":
+        return any(p.type == 'sld' for p in pars)
+    elif kind == "nonmagnetic":
+        return not any(p.type == 'sld' for p in pars)
     return False
 …
     return numpy_dtype, fast, platform
-def list_models_main():
-    # type: () -> None
-    """
-    Run list_models as a main program.  See :func:`list_models` for the
-    kinds of models that can be requested on the command line.
-    """
-    import sys
-    kind = sys.argv[1] if len(sys.argv) > 1 else "all"
-    print("\n".join(list_models(kind)))
 def test_composite_order():
+    """
+    Check that mixture models produce the same result independent of ordder.
+    """
     def test_models(fst, snd):
         """Confirm that two models produce the same parameters"""
 …
     def build_test(first, second):
+        """Construct pair model test"""
         test = lambda description: test_models(first, second)
         description = first + " vs. " + second
 …
     # type: () -> None
     """Check that model load works"""
+    from .product import RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID
     #Test the the model produces the parameters that we would expect
     model = load_model("cylinder@hardsphere*sphere")
     actual = [p.name for p in model.info.parameters.kernel_parameters]
+    target = ("sld sld_solvent radius length theta phi"
+              " radius_effective volfraction "
+              " structure_factor_mode radius_effective_mode"
+              " A_sld A_sld_solvent A_radius").split()
+    target = ["sld", "sld_solvent", "radius", "length", "theta", "phi",
+              RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID,
+              "A_sld", "A_sld_solvent", "A_radius"]
     assert target == actual, "%s != %s"%(target, actual)
+def list_models_main():
+    # type: () -> int
+    """
+    Run list_models as a main program.  See :func:`list_models` for the
+    kinds of models that can be requested on the command line.
+    """
+    import sys
+    kind = sys.argv[1] if len(sys.argv) > 1 else "all"
+    try:
+        models = list_models(kind)
+        print("\n".join(models))
+    except Exception:
+        print(list_models.__doc__)
+        return 1
+    return 0
 if __name__ == "__main__":

sasmodels/data.py

-                      rc88f983
+                      rb297ba9
             mdata = masked_array(data.y, data.mask.copy())
             mdata[~np.isfinite(mdata)] = masked
             if view is 'log':
+            if view == 'log':
                 mdata[mdata <= 0] = masked
             plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.')
 …
             mtheory = masked_array(theory)
             mtheory[~np.isfinite(mtheory)] = masked
             if view is 'log':
+            if view == 'log':
                 mtheory[mtheory <= 0] = masked
             plt.plot(theory_x, scale*mtheory, '-')

sasmodels/details.py

-                      r885753a
+                      rb297ba9
     npars = kernel.info.parameters.npars
     nvalues = kernel.info.parameters.nvalues
     scalars = [value for value, _dispersity, _weight in mesh]
+    scalars = [value for value, dispersity, weight in mesh]
     # skipping scale and background when building values and weights
     _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ())
     #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights)
     length = np.array([len(w) for w in weights])
+    value, dispersity, weight = zip(*mesh[2:npars+2]) if npars else ((), (), ())
+    #weight = correct_theta_weights(kernel.info.parameters, dispersity, weight)
+    length = np.array([len(w) for w in weight])
     offset = np.cumsum(np.hstack((0, length)))
     call_details = make_details(kernel.info, length, offset[:-1], offset[-1])
 …
     data_len = nvalues + 2*sum(len(v) for v in dispersity)
     extra = (32 - data_len%32)%32
     data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra])
+    data = np.hstack((scalars,) + dispersity + weight + ZEROS[:extra])
     data = data.astype(kernel.dtype)
     is_magnetic = convert_magnetism(kernel.info.parameters, data)

sasmodels/direct_model.py

-                      r9150036
+                      re220acc
 from . import resolution2d
 from .details import make_kernel_args, dispersion_mesh
-from .modelinfo import DEFAULT_BACKGROUND
 from .product import RADIUS_MODE_ID
 …
     """
     mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono)
     #print("pars", list(zip(*mesh))[0])
+    #print("in call_kernel: pars:", list(zip(*mesh))[0])
     call_details, values, is_magnetic = make_kernel_args(calculator, mesh)
     #print("values:", values)
+    #print("in call_kernel: values:", values)
     return calculator(call_details, values, cutoff, is_magnetic)
 …
     Use parameter *radius_effective_mode* to select the effective radius
+    calculation.
+    calculation to use amongst the *radius_effective_modes* list given in the
+    model.
     """
     R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0))
     mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono)
     #print("pars", list(zip(*mesh))[0])
+    #print("in call_Fq: pars", list(zip(*mesh))[0])
     call_details, values, is_magnetic = make_kernel_args(calculator, mesh)
     #print("values:", values)
+    #print("in call_Fq: values:", values)
     return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type)
 …
         active = lambda name: True
     #print("pars",[p.id for p in parameters.call_parameters])
+    #print("in get_mesh: pars:",[p.id for p in parameters.call_parameters])
     mesh = [_get_par_weights(p, values, active(p.name))
             for p in parameters.call_parameters]

sasmodels/generate.py

-                      ra8a1d48
+                      ra34b811
     which are hollow.
     *effective_radius(mode, p1, p2, ...)* returns the effective radius of
+    *radius_effective(mode, p1, p2, ...)* returns the effective radius of
     the form with particular dimensions.  Mode determines the type of
     effective radius returned, with mode=1 for equivalent volume.
 …
 def get_data_path(external_dir, target_file):
+    """
+    Search for the target file relative in the installed application.
+    Search first in the location of the generate module in case we are
+    running directly from the distribution.  Search next to the python
+    executable for windows installs.  Search in the ../Resources directory
+    next to the executable for Mac OS/X installs.
+    """
     path = abspath(dirname(__file__))
     if exists(joinpath(path, target_file)):
 …
+#
 # 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.
+# 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|",
 …
 def test_tag_float():
     """check that floating point constants are properly identified and tagged with 'f'"""
+    """Check that floating point constants are identified and tagged with 'f'"""
     cases = """
 …
     return False
+_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE)
+_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]",
+                                   flags=re.MULTILINE)
 def contains_shell_volume(source):
     # type: (List[str]) -> bool
 …
             logger.warn("oriented shapes should define Iqac or Iqabc")
         else:
             raise ValueError("Expected function Iqac or Iqabc for oriented shape")
+            raise ValueError("Expected Iqac or Iqabc for oriented shape")
     # Define the parameter table
 …
                               for p in partable.kernel_parameters))
     # Define the function calls
     call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0"
+    call_radius_effective = "#define CALL_RADIUS_EFFECTIVE(_mode, _v) 0.0"
     if partable.form_volume_parameters:
         refs = _call_pars("_v.", partable.form_volume_parameters)
         if is_hollow:
+            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2)
+            call_volume = (
+                "#define CALL_VOLUME(_form, _shell, _v) "
+                "do { _form = form_volume(%s); _shell = shell_volume(%s); } "
+                "while (0)") % ((",".join(refs),)*2)
         else:
+            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs))
+        if model_info.effective_radius_type:
+            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs))
+            call_volume = (
+                "#define CALL_VOLUME(_form, _shell, _v) "
+                "do { _form = _shell = form_volume(%s); } "
+                "while (0)") % (",".join(refs))
+        if model_info.radius_effective_modes:
+            call_radius_effective = (
+                "#define CALL_RADIUS_EFFECTIVE(_mode, _v) "
+                "radius_effective(_mode, %s)") % (",".join(refs))
     else:
         # Model doesn't have volume.  We could make the kernel run a little
         # faster by not using/transferring the volume normalizations, but
         # the ifdef's reduce readability more than is worthwhile.
+        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)"
+        call_volume = (
+            "#define CALL_VOLUME(_form, _shell, _v) "
+            "do { _form = _shell = 1.0; } while (0)")
     source.append(call_volume)
     source.append(call_effective_radius)
+    source.append(call_radius_effective)
     model_refs = _call_pars("_v.", partable.iq_parameters)
 …
     source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
     source.append("#define PROJECTION %d"%PROJECTION)
+    # TODO: allow mixed python/opencl kernels?
+    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
+    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
+    result = {
+        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
+        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
+    }
+    wrappers = _kernels(kernel_code, call_iq, clear_iq,
+                        call_iqxy, clear_iqxy, model_info.name)
+    code = '\n'.join(source + wrappers[0] + wrappers[1] + wrappers[2])
+    # Note: Identical code for dll and opencl.  This may not always be the case
+    # so leave them as separate entries in the returned value.
+    result = {'dll': code, 'opencl': code}
     return result

sasmodels/gengauss.py

-                      rff31782
+                      rb297ba9
 #!/usr/bin/env python
+"""
+Generate the Gauss-Legendre integration points and save them as a C file.
+"""
 from __future__ import division, print_function
-import sys
 import numpy as np
 …
 def gengauss(n, path):
+    """
+    Save the Gauss-Legendre integration points for length *n* into file *path*.
+    """
     z, w = leggauss(n)

sasmodels/guyou.py

-                      rb3703f5
+                      rb297ba9
 def ellipticJi(u, v, m):
+    """Returns [sn, cn, dn](u + iv|m)."""
     scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m)
     u, v, m = np.broadcast_arrays(u, v, m)
 …
+    ]
-# calculate F(phi+ipsi|m).
-# see Abramowitz and Stegun, 17.4.11.
 def ellipticFi(phi, psi, m):
+    """Returns F(phi+ipsi|m). See Abramowitz and Stegun, 17.4.11."""
     if np.any(phi == 0):
         scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m)
 …
         plt.plot(x, y, 'g')
     for long in range(-limit, limit+1, step):
         x, y = guyou_invert(scale*long, scale*long_line)
+    for longitude in range(-limit, limit+1, step):
+        x, y = guyou_invert(scale*longitude, scale*long_line)
         plt.plot(x, y, 'b')
     plt.xlabel('longitude')

sasmodels/jitter.py

-                      r7d97437
+                      r4e28511
 #!/usr/bin/env python
+# -*- coding: utf-8 -*-
 """
 Jitter Explorer
 …
 Application to explore orientation angle and angular dispersity.
+From the command line::
+    # Show docs
+    python -m sasmodels.jitter --help
+    # Guyou projection jitter, uniform over 20 degree theta and 10 in phi
+    python -m sasmodels.jitter --projection=guyou --dist=uniform --jitter=20,10,0
+From a jupyter cell::
+    import ipyvolume as ipv
+    from sasmodels import jitter
+    import importlib; importlib.reload(jitter)
+    jitter.set_plotter("ipv")
+    size = (10, 40, 100)
+    view = (20, 0, 0)
+    #size = (15, 15, 100)
+    #view = (60, 60, 0)
+    dview = (0, 0, 0)
+    #dview = (5, 5, 0)
+    #dview = (15, 180, 0)
+    #dview = (180, 15, 0)
+    projection = 'equirectangular'
+    #projection = 'azimuthal_equidistance'
+    #projection = 'guyou'
+    #projection = 'sinusoidal'
+    #projection = 'azimuthal_equal_area'
+    dist = 'uniform'
+    #dist = 'gaussian'
+    jitter.run(size=size, view=view, jitter=dview, dist=dist, projection=projection)
+    #filename = projection+('_theta' if dview[0] == 180 else '_phi' if dview[1] == 180 else '')
+    #ipv.savefig(filename+'.png')
 """
 from __future__ import division, print_function
 …
 import argparse
-try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d
-    import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot
-except ImportError:
-    pass
-import matplotlib as mpl
-import matplotlib.pyplot as plt
-from matplotlib.widgets import Slider
 import numpy as np
+from numpy import pi, cos, sin, sqrt, exp, degrees, radians
+def draw_beam(axes, view=(0, 0)):
+from numpy import pi, cos, sin, sqrt, exp, log, degrees, radians, arccos, arctan2
+# Too many complaints about variable names from pylint:
+#    a, b, c, u, v, x, y, z, dx, dy, dz, px, py, pz, R, Rx, Ry, Rz, ...
+# pylint: disable=invalid-name
+def draw_beam(axes, view=(0, 0), alpha=0.5, steps=25):
     """
     Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1)
     """
     #axes.plot([0,0],[0,0],[1,-1])
+    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
+    steps = 25
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(-1, 1, steps)
+    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=alpha)
+    u = np.linspace(0, 2 * pi, steps)
+    v = np.linspace(-1, 1, 2)
     r = 0.02
     x = r*np.outer(np.cos(u), np.ones_like(v))
     y = r*np.outer(np.sin(u), np.ones_like(v))
+    x = r*np.outer(cos(u), np.ones_like(v))
+    y = r*np.outer(sin(u), np.ones_like(v))
     z = 1.3*np.outer(np.ones_like(u), v)
 …
     points = Rz(phi)*Ry(theta)*points
     x, y, z = [v.reshape(shape) for v in points]
+    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
+    axes.plot_surface(x, y, z, color='yellow', alpha=alpha)
+    # TODO: draw endcaps on beam
+    ## Drawing tiny balls on the end will work
+    #draw_sphere(axes, radius=0.02, center=(0, 0, 1.3), color='yellow', alpha=alpha)
+    #draw_sphere(axes, radius=0.02, center=(0, 0, -1.3), color='yellow', alpha=alpha)
+    ## The following does not work
+    #triangles = [(0, i+1, i+2) for i in range(steps-2)]
+    #x_cap, y_cap = x[:, 0], y[:, 0]
+    #for z_cap in z[:, 0], z[:, -1]:
+    #    axes.plot_trisurf(x_cap, y_cap, z_cap, triangles,
+    #                      color='yellow', alpha=alpha)
 def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1):
     """Draw an ellipsoid."""
     a, b, c = size
     u = np.linspace(0, 2 * np.pi, steps)
     v = np.linspace(0, np.pi, steps)
     x = a*np.outer(np.cos(u), np.sin(v))
     y = b*np.outer(np.sin(u), np.sin(v))
     z = c*np.outer(np.ones_like(u), np.cos(v))
+    u = np.linspace(0, 2 * pi, steps)
+    v = np.linspace(0, pi, steps)
+    x = a*np.outer(cos(u), sin(v))
+    y = b*np.outer(sin(u), sin(v))
+    z = c*np.outer(np.ones_like(u), cos(v))
     x, y, z = transform_xyz(view, jitter, x, y, z)
     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
+    axes.plot_surface(x, y, z, color='w', alpha=alpha)
     draw_labels(axes, view, jitter, [
 …
 def draw_sc(axes, size, view, jitter, steps=None, alpha=1):
     """Draw points for simple cubic paracrystal"""
+    # pylint: disable=unused-argument
     atoms = _build_sc()
     _draw_crystal(axes, size, view, jitter, atoms=atoms)
 …
 def draw_fcc(axes, size, view, jitter, steps=None, alpha=1):
     """Draw points for face-centered cubic paracrystal"""
+    # pylint: disable=unused-argument
     # Build the simple cubic crystal
     atoms = _build_sc()
 …
 def draw_bcc(axes, size, view, jitter, steps=None, alpha=1):
     """Draw points for body-centered cubic paracrystal"""
+    # pylint: disable=unused-argument
     # Build the simple cubic crystal
     atoms = _build_sc()
 …
     return atoms
+def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1):
+    """Draw a parallelepiped."""
+def draw_box(axes, size, view):
+    """Draw a wireframe box at a particular view."""
+    a, b, c = size
+    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1])
+    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1])
+    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1])
+    x, y, z = transform_xyz(view, None, x, y, z)
+    def _draw(i, j):
+        axes.plot([x[i], x[j]], [y[i], y[j]], [z[i], z[j]], color='black')
+    _draw(0, 1)
+    _draw(0, 2)
+    _draw(0, 3)
+    _draw(7, 4)
+    _draw(7, 5)
+    _draw(7, 6)
+def draw_parallelepiped(axes, size, view, jitter, steps=None,
+                        color=(0.6, 1.0, 0.6), alpha=1):
+    """Draw a parallelepiped surface, with view and jitter."""
+    # pylint: disable=unused-argument
     a, b, c = size
     x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1])
 …
     x, y, z = transform_xyz(view, jitter, x, y, z)
+    axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
+    # Draw pink face on box.
+    axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha,
+                      linewidth=0)
+    # Colour the c+ face of the box.
     # Since I can't control face color, instead draw a thin box situated just
     # in front of the "c+" face.  Use the c face so that rotations about psi
     # rotate that face.
+    if 1:
+    if 0: # pylint: disable=using-constant-test
+        color = (1, 0.6, 0.6)  # pink
         x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1])
         y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1])
         z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1])
         x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001)
         axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha)
+        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha)
     draw_labels(axes, view, jitter, [
 …
     ])
+def draw_sphere(axes, radius=10., steps=100):
+def draw_sphere(axes, radius=1.0, steps=25,
+                center=(0, 0, 0), color='w', alpha=1.):
     """Draw a sphere"""
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(0, np.pi, steps)
+    x = radius * np.outer(np.cos(u), np.sin(v))
+    y = radius * np.outer(np.sin(u), np.sin(v))
+    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
+    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
+def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0),
+                draw_shape=draw_parallelepiped):
+    u = np.linspace(0, 2 * pi, steps)
+    v = np.linspace(0, pi, steps)
+    x = radius * np.outer(cos(u), sin(v)) + center[0]
+    y = radius * np.outer(sin(u), sin(v)) + center[1]
+    z = radius * np.outer(np.ones(np.size(u)), cos(v)) + center[2]
+    axes.plot_surface(x, y, z, color=color, alpha=alpha)
+    #axes.plot_wireframe(x, y, z)
+def draw_axes(axes, origin=(-1, -1, -1), length=(2, 2, 2)):
+    """Draw wireframe axes lines, with given origin and length"""
+    x, y, z = origin
+    dx, dy, dz = length
+    axes.plot([x, x+dx], [y, y], [z, z], color='black')
+    axes.plot([x, x], [y, y+dy], [z, z], color='black')
+    axes.plot([x, x], [y, y], [z, z+dz], color='black')
+def draw_person_on_sphere(axes, view, height=0.5, radius=1.0):
+    """
+    Draw a person on the surface of a sphere.
+    *view* indicates (latitude, longitude, orientation)
+    """
+    limb_offset = height * 0.05
+    head_radius = height * 0.10
+    head_height = height - head_radius
+    neck_length = head_radius * 0.50
+    shoulder_height = height - 2*head_radius - neck_length
+    torso_length = shoulder_height * 0.55
+    torso_radius = torso_length * 0.30
+    leg_length = shoulder_height - torso_length
+    arm_length = torso_length * 0.90
+    def _draw_part(y, z):
+        x = np.zeros_like(y)
+        xp, yp, zp = transform_xyz(view, None, x, y, z + radius)
+        axes.plot(xp, yp, zp, color='k')
+    # circle for head
+    u = np.linspace(0, 2 * pi, 40)
+    y = head_radius * cos(u)
+    z = head_radius * sin(u) + head_height
+    _draw_part(y, z)
+    # rectangle for body
+    y = np.array([-torso_radius, torso_radius, torso_radius, -torso_radius, -torso_radius])
+    z = np.array([0., 0, torso_length, torso_length, 0]) + leg_length
+    _draw_part(y, z)
+    # arms
+    y = np.array([-torso_radius - limb_offset, -torso_radius - limb_offset, -torso_radius])
+    z = np.array([shoulder_height - arm_length, shoulder_height, shoulder_height])
+    _draw_part(y, z)
+    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type
+    # legs
+    y = np.array([-torso_radius + limb_offset, -torso_radius + limb_offset])
+    z = np.array([0, leg_length])
+    _draw_part(y, z)
+    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type
+    limits = [-radius-height, radius+height]
+    axes.set_xlim(limits)
+    axes.set_ylim(limits)
+    axes.set_zlim(limits)
+    axes.set_axis_off()
+def draw_jitter(axes, view, jitter, dist='gaussian',
+                size=(0.1, 0.4, 1.0),
+                draw_shape=draw_parallelepiped,
+                projection='equirectangular',
+                alpha=0.8,
+                views=None):
     """
     Represent jitter as a set of shapes at different orientations.
     """
+    project, project_weight = get_projection(projection)
     # set max diagonal to 0.95
     scale = 0.95/sqrt(sum(v**2 for v in size))
     size = tuple(scale*v for v in size)
-    #np.random.seed(10)
-    #cloud = np.random.randn(10,3)
-    cloud = [
-        [-1, -1, -1],
-        [-1, -1, +0],
-        [-1, -1, +1],
-        [-1, +0, -1],
-        [-1, +0, +0],
-        [-1, +0, +1],
-        [-1, +1, -1],
-        [-1, +1, +0],
-        [-1, +1, +1],
-        [+0, -1, -1],
-        [+0, -1, +0],
-        [+0, -1, +1],
-        [+0, +0, -1],
-        [+0, +0, +0],
-        [+0, +0, +1],
-        [+0, +1, -1],
-        [+0, +1, +0],
-        [+0, +1, +1],
-        [+1, -1, -1],
-        [+1, -1, +0],
-        [+1, -1, +1],
-        [+1, +0, -1],
-        [+1, +0, +0],
-        [+1, +0, +1],
-        [+1, +1, -1],
-        [+1, +1, +0],
-        [+1, +1, +1],
+    ]
     dtheta, dphi, dpsi = jitter
+    if dtheta == 0:
+        cloud = [v for v in cloud if v[0] == 0]
+    if dphi == 0:
+        cloud = [v for v in cloud if v[1] == 0]
+    if dpsi == 0:
+        cloud = [v for v in cloud if v[2] == 0]
+    draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8)
+    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist]
+    for point in cloud:
+        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]]
+        draw_shape(axes, size, view, delta, alpha=0.8)
+    base = {'gaussian':3, 'rectangle':sqrt(3), 'uniform':1}[dist]
+    def _steps(delta):
+        if views is None:
+            n = max(3, min(25, 2*int(base*delta/5)))
+        else:
+            n = views
+        return base*delta*np.linspace(-1, 1, n) if delta > 0 else [0.]
+    for theta in _steps(dtheta):
+        for phi in _steps(dphi):
+            for psi in _steps(dpsi):
+                w = project_weight(theta, phi, 1.0, 1.0)
+                if w > 0:
+                    dview = project(theta, phi, psi)
+                    draw_shape(axes, size, view, dview, alpha=alpha)
     for v in 'xyz':
         a, b, c = size
         lim = np.sqrt(a**2 + b**2 + c**2)
+        lim = sqrt(a**2 + b**2 + c**2)
         getattr(axes, 'set_'+v+'lim')([-lim, lim])
         getattr(axes, v+'axis').label.set_text(v)
+        #getattr(axes, v+'axis').label.set_text(v)
 PROJECTIONS = [
 …
     'azimuthal_equal_area',
+]
+def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian',
+              projection='equirectangular'):
+    """
+    Draw the dispersion mesh showing the theta-phi orientations at which
+    the model will be evaluated.
+def get_projection(projection):
+    """
     jitter projections
     <https://en.wikipedia.org/wiki/List_of_map_projections>
 …
         Should allow free movement in theta, but phi is distorted.
     """
+    # pylint: disable=unused-argument
     # TODO: try Kent distribution instead of a gaussian warped by projection
-    dist_x = np.linspace(-1, 1, n)
-    weights = np.ones_like(dist_x)
-    if dist == 'gaussian':
-        dist_x *= 3
-        weights = exp(-0.5*dist_x**2)
-    elif dist == 'rectangle':
-        # Note: uses sasmodels ridiculous definition of rectangle width
-        dist_x *= sqrt(3)
-    elif dist == 'uniform':
-        pass
-    else:
-        raise ValueError("expected dist to be gaussian, rectangle or uniform")
     if projection == 'equirectangular':  #define PROJECTION 1
+        def _rotate(theta_i, phi_j):
+            return Rx(phi_j)*Ry(theta_i)
+        def _project(theta_i, phi_j, psi):
+            latitude, longitude = theta_i, phi_j
+            return latitude, longitude, psi, 'xyz'
+            #return Rx(phi_j)*Ry(theta_i)
         def _weight(theta_i, phi_j, w_i, w_j):
             return w_i*w_j*abs(cos(radians(theta_i)))
     elif projection == 'sinusoidal':  #define PROJECTION 2
         def _rotate(theta_i, phi_j):
+        def _project(theta_i, phi_j, psi):
             latitude = theta_i
             scale = cos(radians(latitude))
             longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
+            return Rx(longitude)*Ry(latitude)
+        def _weight(theta_i, phi_j, w_i, w_j):
+            return latitude, longitude, psi, 'xyz'
+            #return Rx(longitude)*Ry(latitude)
+        def _project(theta_i, phi_j, w_i, w_j):
             latitude = theta_i
             scale = cos(radians(latitude))
 …
             return active*w_i*w_j
     elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
         def _rotate(theta_i, phi_j):
+        def _project(theta_i, phi_j, psi):
             from .guyou import guyou_invert
             #latitude, longitude = guyou_invert([theta_i], [phi_j])
             longitude, latitude = guyou_invert([phi_j], [theta_i])
+            return Rx(longitude[0])*Ry(latitude[0])
+            return latitude, longitude, psi, 'xyz'
+            #return Rx(longitude[0])*Ry(latitude[0])
         def _weight(theta_i, phi_j, w_i, w_j):
             return w_i*w_j
+    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
+        def _rotate(theta_i, phi_j):
+    elif projection == 'azimuthal_equidistance':
+        # Note that calculates angles for Rz Ry rather than Rx Ry
+        def _project(theta_i, phi_j, psi):
             latitude = sqrt(theta_i**2 + phi_j**2)
             longitude = degrees(np.arctan2(phi_j, theta_i))
+            longitude = degrees(arctan2(phi_j, theta_i))
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
+            return Rz(longitude)*Ry(latitude)
+            return latitude, longitude, psi-longitude, 'zyz'
+            #R = Rz(longitude)*Ry(latitude)*Rz(psi)
+            #return R_to_xyz(R)
+            #return Rz(longitude)*Ry(latitude)
         def _weight(theta_i, phi_j, w_i, w_j):
             # Weighting for each point comes from the integral:
 …
             return weight*w_i*w_j if latitude < 180 else 0
     elif projection == 'azimuthal_equal_area':
+        def _rotate(theta_i, phi_j):
+        # Note that calculates angles for Rz Ry rather than Rx Ry
+        def _project(theta_i, phi_j, psi):
             radius = min(1, sqrt(theta_i**2 + phi_j**2)/180)
             latitude = 180-degrees(2*np.arccos(radius))
             longitude = degrees(np.arctan2(phi_j, theta_i))
+            latitude = 180-degrees(2*arccos(radius))
+            longitude = degrees(arctan2(phi_j, theta_i))
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
+            return Rz(longitude)*Ry(latitude)
+            return latitude, longitude, psi, 'zyz'
+            #R = Rz(longitude)*Ry(latitude)*Rz(psi)
+            #return R_to_xyz(R)
+            #return Rz(longitude)*Ry(latitude)
         def _weight(theta_i, phi_j, w_i, w_j):
             latitude = sqrt(theta_i**2 + phi_j**2)
 …
         raise ValueError("unknown projection %r"%projection)
+    return _project, _weight
+def R_to_xyz(R):
+    """
+    Return phi, theta, psi Tait-Bryan angles corresponding to the given rotation matrix.
+    Extracting Euler Angles from a Rotation Matrix
+    Mike Day, Insomniac Games
+    https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf
+    Based on: Shoemakeâs "Euler Angle Conversion", Graphics Gems IV, pp.  222-229
+    """
+    phi = arctan2(R[1, 2], R[2, 2])
+    theta = arctan2(-R[0, 2], sqrt(R[0, 0]**2 + R[0, 1]**2))
+    psi = arctan2(R[0, 1], R[0, 0])
+    return degrees(phi), degrees(theta), degrees(psi)
+def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian',
+              projection='equirectangular'):
+    """
+    Draw the dispersion mesh showing the theta-phi orientations at which
+    the model will be evaluated.
+    """
+    _project, _weight = get_projection(projection)
+    def _rotate(theta, phi, z):
+        dview = _project(theta, phi, 0.)
+        if dview[3] == 'zyz':
+            return Rz(dview[1])*Ry(dview[0])*z
+        else:  # dview[3] == 'xyz':
+            return Rx(dview[1])*Ry(dview[0])*z
+    dist_x = np.linspace(-1, 1, n)
+    weights = np.ones_like(dist_x)
+    if dist == 'gaussian':
+        dist_x *= 3
+        weights = exp(-0.5*dist_x**2)
+    elif dist == 'rectangle':
+        # Note: uses sasmodels ridiculous definition of rectangle width
+        dist_x *= sqrt(3)
+    elif dist == 'uniform':
+        pass
+    else:
+        raise ValueError("expected dist to be gaussian, rectangle or uniform")
     # mesh in theta, phi formed by rotating z
     dtheta, dphi, dpsi = jitter
+    dtheta, dphi, dpsi = jitter  # pylint: disable=unused-variable
     z = np.matrix([[0], [0], [radius]])
     points = np.hstack([_rotate(theta_i, phi_j)*z
+    points = np.hstack([_rotate(theta_i, phi_j, z)
                         for theta_i in dtheta*dist_x
                         for phi_j in dphi*dist_x])
 …
     Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
     """
+    dtheta, dphi, dpsi = jitter
+    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    if jitter is None:
+        return points
+    # Hack to deal with the fact that azimuthal_equidistance uses euler angles
+    if len(jitter) == 4:
+        dtheta, dphi, dpsi, _ = jitter
+        points = Rz(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    else:
+        dtheta, dphi, dpsi = jitter
+        points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
     return points
 …
     """
     theta, phi, psi = view
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points # viewing angle
+    #points = Rz(psi)*Ry(pi/2-theta)*Rz(phi)*points # 1-D integration angles
+    #points = Rx(phi)*Ry(theta)*Rz(psi)*points  # angular dispersion angle
     return points
+def orient_relative_to_beam_quaternion(view, points):
+    """
+    Apply the view transform to a set of points.
+    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+    This variant uses quaternions rather than rotation matrices for the
+    computation.  It works but it is not used because it doesn't solve
+    any problems.  The challenge of mapping theta/phi/psi to SO(3) does
+    not disappear by calculating the transform differently.
+    """
+    theta, phi, psi = view
+    x, y, z = [1, 0, 0], [0, 1, 0], [0, 0, 1]
+    q = Quaternion(1, [0, 0, 0])
+    ## Compose a rotation about the three axes by rotating
+    ## the unit vectors before applying the rotation.
+    #q = Quaternion.from_angle_axis(theta, q.rot(x)) * q
+    #q = Quaternion.from_angle_axis(phi, q.rot(y)) * q
+    #q = Quaternion.from_angle_axis(psi, q.rot(z)) * q
+    ## The above turns out to be equivalent to reversing
+    ## the order of application, so ignore it and use below.
+    q = q * Quaternion.from_angle_axis(theta, x)
+    q = q * Quaternion.from_angle_axis(phi, y)
+    q = q * Quaternion.from_angle_axis(psi, z)
+    ## Reverse the order by post-multiply rather than pre-multiply
+    #q = Quaternion.from_angle_axis(theta, x) * q
+    #q = Quaternion.from_angle_axis(phi, y) * q
+    #q = Quaternion.from_angle_axis(psi, z) * q
+    #print("axes psi", q.rot(np.matrix([x, y, z]).T))
+    return q.rot(points)
+#orient_relative_to_beam = orient_relative_to_beam_quaternion
+# === Quaterion class definition === BEGIN
+# Simple stand-alone quaternion class
+# Note: this code works but isn't unused since quaternions didn't solve the
+# representation problem.  Leave it here in case we want to revisit this later.
+#import numpy as np
+class Quaternion(object):
+    r"""
+    Quaternion(w, r) = w + ir[0] + jr[1] + kr[2]
+    Quaternion.from_angle_axis(theta, r) for a rotation of angle theta about
+    an axis oriented toward the direction r.  This defines a unit quaternion,
+    normalizing $r$ to the unit vector $\hat r$, and setting quaternion
+    $Q = \cos \theta + \sin \theta \hat r$
+    Quaternion objects can be multiplied, which applies a rotation about the
+    given axis, allowing composition of rotations without risk of gimbal lock.
+    The resulting quaternion is applied to a set of points using *Q.rot(v)*.
+    """
+    def __init__(self, w, r):
+        self.w = w
+        self.r = np.asarray(r, dtype='d')
+    @staticmethod
+    def from_angle_axis(theta, r):
+        """Build quaternion as rotation theta about axis r"""
+        theta = np.radians(theta)/2
+        r = np.asarray(r)
+        w = np.cos(theta)
+        r = np.sin(theta)*r/np.dot(r, r)
+        return Quaternion(w, r)
+    def __mul__(self, other):
+        """Multiply quaterions"""
+        if isinstance(other, Quaternion):
+            w = self.w*other.w - np.dot(self.r, other.r)
+            r = self.w*other.r + other.w*self.r + np.cross(self.r, other.r)
+            return Quaternion(w, r)
+        raise NotImplementedError("Quaternion * non-quaternion not implemented")
+    def rot(self, v):
+        """Transform point *v* by quaternion"""
+        v = np.asarray(v).T
+        use_transpose = (v.shape[-1] != 3)
+        if use_transpose:
+            v = v.T
+        v = v + np.cross(2*self.r, np.cross(self.r, v) + self.w*v)
+        #v = v + 2*self.w*np.cross(self.r, v) + np.cross(2*self.r, np.cross(self.r, v))
+        if use_transpose:
+            v = v.T
+        return v.T
+    def conj(self):
+        """Conjugate quaternion"""
+        return Quaternion(self.w, -self.r)
+    def inv(self):
+        """Inverse quaternion"""
+        return self.conj()/self.norm()**2
+    def norm(self):
+        """Quaternion length"""
+        return np.sqrt(self.w**2 + np.sum(self.r**2))
+    def __str__(self):
+        return "%g%+gi%+gj%+gk"%(self.w, self.r[0], self.r[1], self.r[2])
+def test_qrot():
+    """Quaternion checks"""
+    # Define rotation of 60 degrees around an axis in y-z that is 60 degrees
+    # from y.  The rotation axis is determined by rotating the point [0, 1, 0]
+    # about x.
+    ax = Quaternion.from_angle_axis(60, [1, 0, 0]).rot([0, 1, 0])
+    q = Quaternion.from_angle_axis(60, ax)
+    # Set the point to be rotated, and its expected rotated position.
+    p = [1, -1, 2]
+    target = [(10+4*np.sqrt(3))/8, (1+2*np.sqrt(3))/8, (14-3*np.sqrt(3))/8]
+    #print(q, q.rot(p) - target)
+    assert max(abs(q.rot(p) - target)) < 1e-14
+#test_qrot()
+#import sys; sys.exit()
+# === Quaterion class definition === END
 # translate between number of dimension of dispersity and the number of
 …
     or the top of the range, depending on whether *mode* is 'central' or 'top'.
     """
+    if portion == 1.0:
+        return data.min(), data.max()
+    elif mode == 'central':
+        data = np.sort(data.flatten())
+        offset = int(portion*len(data)/2 + 0.5)
+        return data[offset], data[-offset]
+    elif mode == 'top':
+        data = np.sort(data.flatten())
+        offset = int(portion*len(data) + 0.5)
+        return data[offset], data[-1]
+    if portion < 1.0:
+        if mode == 'central':
+            data = np.sort(data.flatten())
+            offset = int(portion*len(data)/2 + 0.5)
+            return data[offset], data[-offset]
+        if mode == 'top':
+            data = np.sort(data.flatten())
+            offset = int(portion*len(data) + 0.5)
+            return data[offset], data[-1]
+    # Default: full range
+    return data.min(), data.max()
 def draw_scattering(calculator, axes, view, jitter, dist='gaussian'):
 …
     # compute the pattern
     qx, qy = calculator._data.x_bins, calculator._data.y_bins
+    qx, qy = calculator.qxqy
     Iqxy = calculator(**pars).reshape(len(qx), len(qy))
     # scale it and draw it
     Iqxy = np.log(Iqxy)
+    Iqxy = log(Iqxy)
     if calculator.limits:
         # use limits from orientation (0,0,0)
 …
         vmin = vmax*10**-7
         #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top')
+    #vmin, vmax = Iqxy.min(), Iqxy.max()
     #print("range",(vmin,vmax))
     #qx, qy = np.meshgrid(qx, qy)
     if 0:
+    if 0:  # pylint: disable=using-constant-test
         level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i')
         level[level < 0] = 0
+        from matplotlib import pylab as plt
         colors = plt.get_cmap()(level)
+        axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
+    elif 1:
+        #from matplotlib import cm
+        #colors = cm.coolwarm(level)
+        #colors = cm.gist_yarg(level)
+        #colors = cm.Wistia(level)
+        colors[level <= 0, 3] = 0.  # set floor to transparent
+        x, y = np.meshgrid(qx/qx.max(), qy/qy.max())
+        axes.plot_surface(x, y, -1.1*np.ones_like(x), facecolors=colors)
+    elif 1:  # pylint: disable=using-constant-test
         axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
                       levels=np.linspace(vmin, vmax, 24))
 …
     calculator = DirectModel(data, model)
+    # Remember the data axes so we can plot the results
+    calculator.qxqy = (q, q)
     # stuff the values for non-orientation parameters into the calculator
     calculator.pars = pars.copy()
 …
     # under rotation or angular dispersion
     Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars)
     Iqxy = np.log(Iqxy)
+    Iqxy = log(Iqxy)
     vmin, vmax = clipped_range(Iqxy, 0.95, mode='top')
     calculator.limits = vmin, vmax+1
 …
+}
 def run(model_name='parallelepiped', size=(10, 40, 100),
+        view=(0, 0, 0), jitter=(0, 0, 0),
         dist='gaussian', mesh=30,
         projection='equirectangular'):
 …
     *size* gives the dimensions (a, b, c) of the shape.
+    *view* gives the initial view (theta, phi, psi) of the shape.
+    *view* gives the initial jitter (dtheta, dphi, dpsi) of the shape.
     *dist* is the type of dispersition: gaussian, rectangle, or uniform.
 …
     calculator, size = select_calculator(model_name, n=150, size=size)
     draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped)
+    #draw_shape = draw_fcc
     ## uncomment to set an independent the colour range for every view
 …
     calculator.limits = None
+    ## initial view
+    #theta, dtheta = 70., 10.
+    #phi, dphi = -45., 3.
+    #psi, dpsi = -45., 3.
+    theta, phi, psi = 0, 0, 0
+    dtheta, dphi, dpsi = 0, 0, 0
+    PLOT_ENGINE(calculator, draw_shape, size, view, jitter, dist, mesh, projection)
+def _mpl_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection):
+    # Note: travis-ci does not support mpl_toolkits.mplot3d, but this shouldn't be
+    # an issue since we are lazy-loading the package on a path that isn't tested.
+    # Importing mplot3d adds projection='3d' option to subplot
+    import mpl_toolkits.mplot3d  # pylint: disable=unused-variable
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    from matplotlib.widgets import Slider
     ## create the plot window
 …
     ## add control widgets to plot
+    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props)
+    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props)
+    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props)
+    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta)
+    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi)
+    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi)
+    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props)
+    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props)
+    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props)
+    axes_theta = plt.axes([0.05, 0.15, 0.50, 0.04], **props)
+    axes_phi = plt.axes([0.05, 0.10, 0.50, 0.04], **props)
+    axes_psi = plt.axes([0.05, 0.05, 0.50, 0.04], **props)
+    stheta = Slider(axes_theta, u'Îž', -90, 90, valinit=0)
+    sphi = Slider(axes_phi, u'Ï', -180, 180, valinit=0)
+    spsi = Slider(axes_psi, u'Ï', -180, 180, valinit=0)
+    axes_dtheta = plt.axes([0.70, 0.15, 0.20, 0.04], **props)
+    axes_dphi = plt.axes([0.70, 0.1, 0.20, 0.04], **props)
+    axes_dpsi = plt.axes([0.70, 0.05, 0.20, 0.04], **props)
     # Note: using ridiculous definition of rectangle distribution, whose width
     # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep
     # the maximum width to 90.
     dlimit = DIST_LIMITS[dist]
+    sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta)
+    sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi)
+    sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi)
+    sdtheta = Slider(axes_dtheta, u'ÎÎž', 0, 2*dlimit, valinit=0)
+    sdphi = Slider(axes_dphi, u'ÎÏ', 0, 2*dlimit, valinit=0)
+    sdpsi = Slider(axes_dpsi, u'ÎÏ', 0, 2*dlimit, valinit=0)
+    ## initial view and jitter
+    theta, phi, psi = view
+    stheta.set_val(theta)
+    sphi.set_val(phi)
+    spsi.set_val(psi)
+    dtheta, dphi, dpsi = jitter
+    sdtheta.set_val(dtheta)
+    sdphi.set_val(dphi)
+    sdpsi.set_val(dpsi)
     ## callback to draw the new view
+    def update(val, axis=None):
+    def _update(val, axis=None):
+        # pylint: disable=unused-argument
         view = stheta.val, sphi.val, spsi.val
         jitter = sdtheta.val, sdphi.val, sdpsi.val
         # set small jitter as 0 if multiple pd dims
         dims = sum(v > 0 for v in jitter)
         limit = [0, 0.5, 5][dims]
+        limit = [0, 0.5, 5, 5][dims]
         jitter = [0 if v < limit else v for v in jitter]
         axes.cla()
+        draw_beam(axes, (0, 0))
+        draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape)
+        #draw_jitter(axes, view, (0,0,0))
+        ## Visualize as person on globe
+        #draw_sphere(axes, radius=0.5)
+        #draw_person_on_sphere(axes, view, radius=0.5)
+        ## Move beam instead of shape
+        #draw_beam(axes, -view[:2])
+        #draw_jitter(axes, (0,0,0), (0,0,0), views=3)
+        ## Move shape and draw scattering
+        draw_beam(axes, (0, 0), alpha=1.)
+        #draw_person_on_sphere(axes, view, radius=1.2, height=0.5)
+        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.,
+                    draw_shape=draw_shape, projection=projection, views=3)
         draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection)
         draw_scattering(calculator, axes, view, jitter, dist=dist)
         plt.gcf().canvas.draw()
     ## bind control widgets to view updater
     stheta.on_changed(lambda v: update(v, 'theta'))
     sphi.on_changed(lambda v: update(v, 'phi'))
     spsi.on_changed(lambda v: update(v, 'psi'))
     sdtheta.on_changed(lambda v: update(v, 'dtheta'))
     sdphi.on_changed(lambda v: update(v, 'dphi'))
     sdpsi.on_changed(lambda v: update(v, 'dpsi'))
+    stheta.on_changed(lambda v: _update(v, 'theta'))
+    sphi.on_changed(lambda v: _update(v, 'phi'))
+    spsi.on_changed(lambda v: _update(v, 'psi'))
+    sdtheta.on_changed(lambda v: _update(v, 'dtheta'))
+    sdphi.on_changed(lambda v: _update(v, 'dphi'))
+    sdpsi.on_changed(lambda v: _update(v, 'dpsi'))
     ## initialize view
     update(None, 'phi')
+    _update(None, 'phi')
     ## go interactive
     plt.show()
+def map_colors(z, kw):
+    """
+    Process matplotlib-style colour arguments.
+    Pulls 'cmap', 'alpha', 'vmin', and 'vmax' from th *kw* dictionary, setting
+    the *kw['color']* to an RGB array.  These are ignored if 'c' or 'color' are
+    set inside *kw*.
+    """
+    from matplotlib import cm
+    cmap = kw.pop('cmap', cm.coolwarm)
+    alpha = kw.pop('alpha', None)
+    vmin = kw.pop('vmin', z.min())
+    vmax = kw.pop('vmax', z.max())
+    c = kw.pop('c', None)
+    color = kw.pop('color', c)
+    if color is None:
+        znorm = ((z - vmin) / (vmax - vmin)).clip(0, 1)
+        color = cmap(znorm)
+    elif isinstance(color, np.ndarray) and color.shape == z.shape:
+        color = cmap(color)
+    if alpha is None:
+        if isinstance(color, np.ndarray):
+            color = color[..., :3]
+    else:
+        color[..., 3] = alpha
+    kw['color'] = color
+def make_vec(*args):
+    """Turn all elements of *args* into numpy arrays"""
+    #return [np.asarray(v, 'd').flatten() for v in args]
+    return [np.asarray(v, 'd') for v in args]
+def make_image(z, kw):
+    """Convert numpy array *z* into a *PIL* RGB image."""
+    import PIL.Image
+    from matplotlib import cm
+    cmap = kw.pop('cmap', cm.coolwarm)
+    znorm = (z-z.min())/z.ptp()
+    c = cmap(znorm)
+    c = c[..., :3]
+    rgb = np.asarray(c*255, 'u1')
+    image = PIL.Image.fromarray(rgb, mode='RGB')
+    return image
+_IPV_MARKERS = {
+    'o': 'sphere',
+}
+_IPV_COLORS = {
+    'w': 'white',
+    'k': 'black',
+    'c': 'cyan',
+    'm': 'magenta',
+    'y': 'yellow',
+    'r': 'red',
+    'g': 'green',
+    'b': 'blue',
+}
+def _ipv_fix_color(kw):
+    alpha = kw.pop('alpha', None)
+    color = kw.get('color', None)
+    if isinstance(color, str):
+        color = _IPV_COLORS.get(color, color)
+        kw['color'] = color
+    if alpha is not None:
+        color = kw['color']
+        #TODO: convert color to [r, g, b, a] if not already
+        if isinstance(color, (tuple, list)):
+            if len(color) == 3:
+                color = (color[0], color[1], color[2], alpha)
+            else:
+                color = (color[0], color[1], color[2], alpha*color[3])
+            color = np.array(color)
+        if isinstance(color, np.ndarray) and color.shape[-1] == 4:
+            color[..., 3] = alpha
+            kw['color'] = color
+def _ipv_set_transparency(kw, obj):
+    color = kw.get('color', None)
+    if (isinstance(color, np.ndarray)
+            and color.shape[-1] == 4
+            and (color[..., 3] != 1.0).any()):
+        obj.material.transparent = True
+        obj.material.side = "FrontSide"
+def ipv_axes():
+    """
+    Build a matplotlib style Axes interface for ipyvolume
+    """
+    import ipyvolume as ipv
+    class Axes(object):
+        """
+        Matplotlib Axes3D style interface to ipyvolume renderer.
+        """
+        # pylint: disable=no-self-use,no-init
+        # transparency can be achieved by setting the following:
+        #    mesh.color = [r, g, b, alpha]
+        #    mesh.material.transparent = True
+        #    mesh.material.side = "FrontSide"
+        # smooth(ish) rotation can be achieved by setting:
+        #    slide.continuous_update = True
+        #    figure.animation = 0.
+        #    mesh.material.x = x
+        #    mesh.material.y = y
+        #    mesh.material.z = z
+        # maybe need to synchronize update of x/y/z to avoid shimmy when moving
+        def plot(self, x, y, z, **kw):
+            """mpl style plot interface for ipyvolume"""
+            _ipv_fix_color(kw)
+            x, y, z = make_vec(x, y, z)
+            ipv.plot(x, y, z, **kw)
+        def plot_surface(self, x, y, z, **kw):
+            """mpl style plot_surface interface for ipyvolume"""
+            facecolors = kw.pop('facecolors', None)
+            if facecolors is not None:
+                kw['color'] = facecolors
+            _ipv_fix_color(kw)
+            x, y, z = make_vec(x, y, z)
+            h = ipv.plot_surface(x, y, z, **kw)
+            _ipv_set_transparency(kw, h)
+            #h.material.side = "DoubleSide"
+            return h
+        def plot_trisurf(self, x, y, triangles=None, Z=None, **kw):
+            """mpl style plot_trisurf interface for ipyvolume"""
+            kw.pop('linewidth', None)
+            _ipv_fix_color(kw)
+            x, y, z = make_vec(x, y, Z)
+            if triangles is not None:
+                triangles = np.asarray(triangles)
+            h = ipv.plot_trisurf(x, y, z, triangles=triangles, **kw)
+            _ipv_set_transparency(kw, h)
+            return h
+        def scatter(self, x, y, z, **kw):
+            """mpl style scatter interface for ipyvolume"""
+            x, y, z = make_vec(x, y, z)
+            map_colors(z, kw)
+            marker = kw.get('marker', None)
+            kw['marker'] = _IPV_MARKERS.get(marker, marker)
+            h = ipv.scatter(x, y, z, **kw)
+            _ipv_set_transparency(kw, h)
+            return h
+        def contourf(self, x, y, v, zdir='z', offset=0, levels=None, **kw):
+            """mpl style contourf interface for ipyvolume"""
+            # pylint: disable=unused-argument
+            # Don't use contour for now (although we might want to later)
+            self.pcolor(x, y, v, zdir='z', offset=offset, **kw)
+        def pcolor(self, x, y, v, zdir='z', offset=0, **kw):
+            """mpl style pcolor interface for ipyvolume"""
+            # pylint: disable=unused-argument
+            x, y, v = make_vec(x, y, v)
+            image = make_image(v, kw)
+            xmin, xmax = x.min(), x.max()
+            ymin, ymax = y.min(), y.max()
+            x = np.array([[xmin, xmax], [xmin, xmax]])
+            y = np.array([[ymin, ymin], [ymax, ymax]])
+            z = x*0 + offset
+            u = np.array([[0., 1], [0, 1]])
+            v = np.array([[0., 0], [1, 1]])
+            h = ipv.plot_mesh(x, y, z, u=u, v=v, texture=image, wireframe=False)
+            _ipv_set_transparency(kw, h)
+            h.material.side = "DoubleSide"
+            return h
+        def text(self, *args, **kw):
+            """mpl style text interface for ipyvolume"""
+            pass
+        def set_xlim(self, limits):
+            """mpl style set_xlim interface for ipyvolume"""
+            ipv.xlim(*limits)
+        def set_ylim(self, limits):
+            """mpl style set_ylim interface for ipyvolume"""
+            ipv.ylim(*limits)
+        def set_zlim(self, limits):
+            """mpl style set_zlim interface for ipyvolume"""
+            ipv.zlim(*limits)
+        def set_axes_on(self):
+            """mpl style set_axes_on interface for ipyvolume"""
+            ipv.style.axis_on()
+        def set_axis_off(self):
+            """mpl style set_axes_off interface for ipyvolume"""
+            ipv.style.axes_off()
+    return Axes()
+def _ipv_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection):
+    from IPython.display import display
+    import ipywidgets as widgets
+    import ipyvolume as ipv
+    axes = ipv_axes()
+    def _draw(view, jitter):
+        camera = ipv.gcf().camera
+        #print(ipv.gcf().__dict__.keys())
+        #print(dir(ipv.gcf()))
+        ipv.figure(animation=0.)  # no animation when updating object mesh
+        # set small jitter as 0 if multiple pd dims
+        dims = sum(v > 0 for v in jitter)
+        limit = [0, 0.5, 5, 5][dims]
+        jitter = [0 if v < limit else v for v in jitter]
+        ## Visualize as person on globe
+        #draw_beam(axes, (0, 0))
+        #draw_sphere(axes, radius=0.5)
+        #draw_person_on_sphere(axes, view, radius=0.5)
+        ## Move beam instead of shape
+        #draw_beam(axes, view=(-view[0], -view[1]))
+        #draw_jitter(axes, view=(0,0,0), jitter=(0,0,0))
+        ## Move shape and draw scattering
+        draw_beam(axes, (0, 0), steps=25)
+        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.0,
+                    draw_shape=draw_shape, projection=projection)
+        draw_mesh(axes, view, jitter, dist=dist, n=mesh, radius=0.95,
+                  projection=projection)
+        draw_scattering(calculator, axes, view, jitter, dist=dist)
+        draw_axes(axes, origin=(-1, -1, -1.1))
+        ipv.style.box_off()
+        ipv.style.axes_off()
+        ipv.xyzlabel(" ", " ", " ")
+        ipv.gcf().camera = camera
+        ipv.show()
+    trange, prange = (-180., 180., 1.), (-180., 180., 1.)
+    dtrange, dprange = (0., 180., 1.), (0., 180., 1.)
+    ## Super simple interfaca, but uses non-ascii variable namese
+    # Îž Ï Ï ÎÎž ÎÏ ÎÏ
+    #def update(**kw):
+    #    view = kw['Îž'], kw['Ï'], kw['Ï']
+    #    jitter = kw['ÎÎž'], kw['ÎÏ'], kw['ÎÏ']
+    #    draw(view, jitter)
+    #widgets.interact(update, Îž=trange, Ï=prange, Ï=prange, ÎÎž=dtrange, ÎÏ=dprange, ÎÏ=dprange)
+    def _update(theta, phi, psi, dtheta, dphi, dpsi):
+        _draw(view=(theta, phi, psi), jitter=(dtheta, dphi, dpsi))
+    def _slider(name, slice, init=0.):
+        return widgets.FloatSlider(
+            value=init,
+            min=slice[0],
+            max=slice[1],
+            step=slice[2],
+            description=name,
+            disabled=False,
+            #continuous_update=True,
+            continuous_update=False,
+            orientation='horizontal',
+            readout=True,
+            readout_format='.1f',
+            )
+    theta = _slider(u'Îž', trange, view[0])
+    phi = _slider(u'Ï', prange, view[1])
+    psi = _slider(u'Ï', prange, view[2])
+    dtheta = _slider(u'ÎÎž', dtrange, jitter[0])
+    dphi = _slider(u'ÎÏ', dprange, jitter[1])
+    dpsi = _slider(u'ÎÏ', dprange, jitter[2])
+    fields = {
+        'theta': theta, 'phi': phi, 'psi': psi,
+        'dtheta': dtheta, 'dphi': dphi, 'dpsi': dpsi,
+    }
+    ui = widgets.HBox([
+        widgets.VBox([theta, phi, psi]),
+        widgets.VBox([dtheta, dphi, dpsi])
+    ])
+    out = widgets.interactive_output(_update, fields)
+    display(ui, out)
+_ENGINES = {
+    "matplotlib": _mpl_plot,
+    "mpl": _mpl_plot,
+    #"plotly": _plotly_plot,
+    "ipvolume": _ipv_plot,
+    "ipv": _ipv_plot,
+}
+PLOT_ENGINE = _ENGINES["matplotlib"]
+def set_plotter(name):
+    """
+    Setting the plotting engine to matplotlib/ipyvolume or equivalently mpl/ipv.
+    """
+    global PLOT_ENGINE
+    PLOT_ENGINE = _ENGINES[name]
 def main():
+    """
+    Command line interface to the jitter viewer.
+    """
     parser = argparse.ArgumentParser(
         description="Display jitter",
 …
     parser.add_argument('-s', '--size', type=str, default='10,40,100',
                         help='a,b,c lengths')
+    parser.add_argument('-v', '--view', type=str, default='0,0,0',
+                        help='initial view angles')
+    parser.add_argument('-j', '--jitter', type=str, default='0,0,0',
+                        help='initial angular dispersion')
     parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS,
                         default=DISTRIBUTIONS[0],
 …
                         help='oriented shape')
     opts = parser.parse_args()
+    size = tuple(int(v) for v in opts.size.split(','))
+    run(opts.shape, size=size,
+    size = tuple(float(v) for v in opts.size.split(','))
+    view = tuple(float(v) for v in opts.view.split(','))
+    jitter = tuple(float(v) for v in opts.jitter.split(','))
+    run(opts.shape, size=size, view=view, jitter=jitter,
         mesh=opts.mesh, dist=opts.distribution,
         projection=opts.projection)

sasmodels/kernel.py

-                      rcd28947
+                      ra34b811
 # pylint: enable=unused-import
 class KernelModel(object):
+    """
+    Model definition for the compute engine.
+    """
     info = None  # type: ModelInfo
     dtype = None # type: np.dtype
     def make_kernel(self, q_vectors):
         # type: (List[np.ndarray]) -> "Kernel"
+        """
+        Instantiate a kernel for evaluating the model at *q_vectors*.
+        """
         raise NotImplementedError("need to implement make_kernel")
     def release(self):
         # type: () -> None
+        """
+        Free resources associated with the kernel.
+        """
         pass
 class Kernel(object):
+    #: kernel dimension, either "1d" or "2d"
+    """
+    Instantiated model for the compute engine, applied to a particular *q*.
+    Subclasses should define :meth:`_call_kernel` to evaluate the kernel over
+    its inputs.
+    """
+    #: Kernel dimension, either "1d" or "2d".
     dim = None  # type: str
+    #: Model info.
     info = None  # type: ModelInfo
     results = None # type: List[np.ndarray]
+    #: Numerical precision for the computation.
     dtype = None  # type: np.dtype
+    #: q values at which the kernel is to be evaluated
+    q_input = None  # type: Any
+    #: Place to hold result of :meth:`_call_kernel` for subclass.
+    result = None # type: np.ndarray
     def Iq(self, call_details, values, cutoff, magnetic):
 …
         built into the model, and do not need an additional scale.
         """
         _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic,
                               effective_radius_type=0)
+        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff,
+                                            magnetic, radius_effective_mode=0)
         combined_scale = values[0]/shell_volume
         background = values[1]
 …
     __call__ = Iq
+    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0):
+    def Fq(self, call_details, values, cutoff, magnetic,
+           radius_effective_mode=0):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
         r"""
 …
         .. math::
             I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background}
+            I(q) = \text{scale} P (1 + <F>^2/<F^2> (S - 1)) + \text{background}
                  = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background}
 …
         .. math::
             P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k}
+            P(q)=\frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k}
         The form factor itself is scaled by volume and contrast to compute the
 …
         volume fraction of the particles.  The model can have several
         different ways to compute effective radius, with the
         *effective_radius_type* parameter used to select amongst them.  The
+        *radius_effective_mode* parameter used to select amongst them.  The
         volume fraction of particles should be determined from the total
         volume fraction of the form, not just the shell volume fraction.
 …
         hollow and solid shapes.
         """
+        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type)
+        self._call_kernel(call_details, values, cutoff, magnetic,
+                          radius_effective_mode)
         #print("returned",self.q_input.q, self.result)
         nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
 …
         form_volume = self.result[nout*self.q_input.nq + 1]/total_weight
         shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight
         effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight
+        radius_effective = self.result[nout*self.q_input.nq + 3]/total_weight
         if shell_volume == 0.:
             shell_volume = 1.
+        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None
+        F1 = (self.result[1:nout*self.q_input.nq:nout]/total_weight
+              if nout == 2 else None)
         F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight
         return F1, F2, effective_radius, shell_volume, form_volume/shell_volume
+        return F1, F2, radius_effective, shell_volume, form_volume/shell_volume
     def release(self):
         # type: () -> None
+        """
+        Free resources associated with the kernel instance.
+        """
         pass
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+    def _call_kernel(self, call_details, values, cutoff, magnetic,
+                     radius_effective_mode):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
         """

sasmodels/kernel_iq.c

-                      r12f4c19
+                      rde032da
 //      parameters in the parameter table.
 //  CALL_VOLUME(form, shell, table) : assign form and shell values
 //  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function
+//  CALL_RADIUS_EFFECTIVE(mode, table) : call the R_eff function
 //  CALL_IQ(q, table) : call the Iq function for 1D calcs.
 //  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data.
 …
 static void set_spin_weights(double in_spin, double out_spin, double weight[6])
+{
+  double norm;
   in_spin = clip(in_spin, 0.0, 1.0);
   out_spin = clip(out_spin, 0.0, 1.0);
 …
   // However, since the weights are applied to the final intensity and
   // are not interned inside the I(q) function, we want the full
+  // weight and not the square root.  Any function using
+  // set_spin_weights as part of calculating an amplitude will need to
+  // manually take that square root, but there is currently no such
+  // function.
+  weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd
+  weight[1] = (1.0-in_spin) * out_spin;       // du
+  weight[2] = in_spin * (1.0-out_spin);       // ud
+  weight[3] = in_spin * out_spin;             // uu
+  // weight and not the square root.  Anyway no function will ever use
+  // set_spin_weights as part of calculating an amplitude, as the weights are
+  // related to polarisation efficiency of the instrument. The weights serve to
+  // construct various magnet scattering cross sections, which are linear combinations
+  // of the spin-resolved cross sections. The polarisation efficiency e_in and e_out
+  // are parameters ranging from 0.5 (unpolarised) beam to 1 (perfect optics).
+  // For in_spin or out_spin <0.5 one assumes a CS, where the spin is reversed/flipped
+  // with respect to the initial supermirror polariser. The actual polarisation efficiency
+  // in this case is however e_in/out = 1-in/out_spin.
+  if (out_spin < 0.5){norm=1-out_spin;}
+  else{norm=out_spin;}
+// The norm is needed to make sure that the scattering cross sections are
+//correctly weighted, such that the sum of spin-resolved measurements adds up to
+// the unpolarised or half-polarised scattering cross section. No intensity weighting
+// needed on the incoming polariser side (assuming that a user), has normalised
+// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively.
+  weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd
+  weight[1] = (1.0-in_spin) * out_spin / norm;       // du
+  weight[2] = in_spin * (1.0-out_spin) / norm;       // ud
+  weight[3] = in_spin * out_spin / norm;             // uu
   weight[4] = weight[1]; // du.imag
   weight[5] = weight[2]; // ud.imag
 …
     switch (xs) {
       default: // keep compiler happy; condition ensures xs in [0,1,2,3]
       case 0: // uu => sld - D M_perpx
+      case 0: // dd => sld - D M_perpx
           return sld - px*perp;
       case 1: // ud.real => -D M_perpy
+      case 1: // du.real => -D M_perpy
           return py*perp;
       case 2: // du.real => -D M_perpy
+      case 2: // ud.real => -D M_perpy
           return py*perp;
       case 3: // dd => sld + D M_perpx
+      case 3: // uu => sld + D M_perpx
           return sld + px*perp;
+    }
 …
     pglobal double *result,       // nq+1 return values, again with padding
     const double cutoff,          // cutoff in the dispersity weight product
     int32_t effective_radius_type // which effective radius to compute
+    int32_t radius_effective_mode // which effective radius to compute
+    )
+{
 …
       weighted_form += weight * form;
       weighted_shell += weight * shell;
       if (effective_radius_type != 0) {
         weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table);
+      if (radius_effective_mode != 0) {
+        weighted_radius += weight * CALL_RADIUS_EFFECTIVE(radius_effective_mode, local_values.table);
+      }
       BUILD_ROTATION();

sasmodels/kernelcl.py

-                      r8064d5e
+                      ra34b811
 from __future__ import print_function
+import sys
 import os
 import warnings
 …
     from time import perf_counter as clock
 except ImportError: # CRUFT: python < 3.3
-    import sys
     if sys.platform.count("darwin") > 0:
         from time import time as clock
 …
 import numpy as np  # type: ignore
 # Attempt to setup opencl. This may fail if the pyopencl package is not
+# Attempt to setup OpenCL. This may fail if the pyopencl package is not
 # installed or if it is installed but there are no devices available.
 try:
 …
     from pyopencl import mem_flags as mf
     from pyopencl.characterize import get_fast_inaccurate_build_options
     # Ask OpenCL for the default context so that we know that one exists
+    # Ask OpenCL for the default context so that we know that one exists.
     cl.create_some_context(interactive=False)
     HAVE_OPENCL = True
 …
 # pylint: enable=unused-import
+# CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path)
+# CRUFT: pyopencl < 2017.1 (as of June 2016 needs quotes around include path).
 def quote_path(v):
+    # type: (str) -> str
     """
     Quote the path if it is not already quoted.
 …
     return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v
 def fix_pyopencl_include():
+    # type: (None) -> None
     """
     Monkey patch pyopencl to allow spaces in include file path.
     """
+    import pyopencl as cl
+    if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'):
+        cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS]
+    # pylint: disable=protected-access
+    import pyopencl
+    if hasattr(pyopencl, '_DEFAULT_INCLUDE_OPTIONS'):
+        pyopencl._DEFAULT_INCLUDE_OPTIONS = [
+            quote_path(v) for v in pyopencl._DEFAULT_INCLUDE_OPTIONS
+            ]
 if HAVE_OPENCL:
 …
 MAX_LOOPS = 2048
 # Pragmas for enable OpenCL features.  Be sure to protect them so that they
 # still compile even if OpenCL is not present.
 …
 """
 def use_opencl():
+    # type: () -> bool
+    """Return True if OpenCL is the default computational engine"""
     sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower()
     return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda")
 ENV = None
 def reset_environment():
+    # type: () -> None
     """
     Call to create a new OpenCL context, such as after a change to SAS_OPENCL.
 …
     global ENV
     ENV = GpuEnvironment() if use_opencl() else None
 def environment():
 …
     return ENV
 def has_type(device, dtype):
     # type: (cl.Device, np.dtype) -> bool
 …
         return "cl_khr_fp64" in device.extensions
     else:
         # Not supporting F16 type since it isn't accurate enough
+        # Not supporting F16 type since it isn't accurate enough.
         return False
 def get_warp(kernel, queue):
 …
         cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
         queue.device)
 def compile_model(context, source, dtype, fast=False):
 …
         source_list.insert(0, _F64_PRAGMA)
     # Note: USE_SINCOS makes the intel cpu slower under opencl
+    # Note: USE_SINCOS makes the Intel CPU slower under OpenCL.
     if context.devices[0].type == cl.device_type.GPU:
         source_list.insert(0, "#define USE_SINCOS\n")
 …
     source = "\n".join(source_list)
     program = cl.Program(context, source).build(options=options)
     #print("done with "+program)
     return program
 # for now, this returns one device in the context
 # TODO: create a context that contains all devices on all platforms
+# For now, this returns one device in the context.
+# TODO: Create a context that contains all devices on all platforms.
 class GpuEnvironment(object):
     """
+    GPU context, with possibly many devices, and one queue per device.
+    Because the environment can be reset during a live program (e.g., if the
+    user changes the active GPU device in the GUI), everything associated
+    with the device context must be cached in the environment and recreated
+    if the environment changes.  The *cache* attribute is a simple dictionary
+    which holds keys and references to objects, such as compiled kernels and
+    allocated buffers.  The running program should check in the cache for
+    long lived objects and create them if they are not there.  The program
+    should not hold onto cached objects, but instead only keep them active
+    for the duration of a function call.  When the environment is destroyed
+    then the *release* method for each active cache item is called before
+    the environment is freed.  This means that each cl buffer should be
+    in its own cache entry.
+    GPU context for OpenCL, with possibly many devices and one queue per device.
     """
     def __init__(self):
         # type: () -> None
         # find gpu context
+        # Find gpu context.
         context_list = _create_some_context()
 …
                 self.context[dtype] = None
         # Build a queue for each context
+        # Build a queue for each context.
         self.queue = {}
         context = self.context[F32]
 …
             self.queue[F64] = cl.CommandQueue(context, context.devices[0])
         # Byte boundary for data alignment
+        ## Byte boundary for data alignment.
         #self.data_boundary = max(context.devices[0].min_data_type_align_size
         #                         for context in self.context.values())
         # Cache for compiled programs, and for items in context
+        # Cache for compiled programs, and for items in context.
         self.compiled = {}
-        self.cache = {}
     def has_type(self, dtype):
 …
         """
         # Note: PyOpenCL caches based on md5 hash of source, options and device
+        # so we don't really need to cache things for ourselves.  I'll do so
+        # anyway just to save some data munging time.
+        # but I'll do so as well just to save some data munging time.
         tag = generate.tag_source(source)
         key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else ""))
         # Check timestamp on program
+        # Check timestamp on program.
         program, program_timestamp = self.compiled.get(key, (None, np.inf))
         if program_timestamp < timestamp:
 …
         return program
-    def free_buffer(self, key):
-        if key in self.cache:
-            self.cache[key].release()
-            del self.cache[key]
-    def __del__(self):
-        for v in self.cache.values():
-            release = getattr(v, 'release', lambda: None)
-            release()
-        self.cache = {}
-_CURRENT_ID = 0
-def unique_id():
-    global _CURRENT_ID
-    _CURRENT_ID += 1
-    return _CURRENT_ID
 def _create_some_context():
 …
     which one (and not a CUDA device, or no GPU).
     """
     # Assume we do not get here if SAS_OPENCL is None or CUDA
+    # Assume we do not get here if SAS_OPENCL is None or CUDA.
     sas_opencl = os.environ.get('SAS_OPENCL', 'opencl')
     if sas_opencl.lower() != 'opencl':
         # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context
+        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context.
         os.environ["PYOPENCL_CTX"] = sas_opencl
 …
             return [cl.create_some_context(interactive=False)]
         except Exception as exc:
+            # TODO: Should warnings instead be put into logging.warn?
             warnings.warn(str(exc))
+            warnings.warn("pyopencl.create_some_context() failed")
+            warnings.warn("the environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set correctly")
+            warnings.warn(
+                "pyopencl.create_some_context() failed.  The environment "
+                "variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set "
+                "correctly")
     return _get_default_context()
 def _get_default_context():
 …
     # is running may increase throughput.
+    #
     # Macbook pro, base install:
+    # MacBook Pro, base install:
     #     {'Apple': [Intel CPU, NVIDIA GPU]}
     # Macbook pro, base install:
+    # MacBook Pro, base install:
     #     {'Apple': [Intel CPU, Intel GPU]}
     # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed
+    # 2 x NVIDIA 295 with Intel and NVIDIA opencl drivers install:
     #     {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]}
     gpu, cpu = None, None
 …
             else:
                 # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM
                 # Intel Phi for example registers as an accelerator
+                # Intel Phi for example registers as an accelerator.
                 # Since the user installed a custom device on their system
                 # and went through the pain of sorting out OpenCL drivers for
 …
                 gpu = device
     # order the devices by gpu then by cpu; when searching for an available
+    # Order the devices by gpu then by cpu; when searching for an available
     # device by data type they will be checked in this order, which means
     # that if the gpu supports double then the cpu will never be used (though
 …
     that the compiler is allowed to take shortcuts.
     """
+    info = None  # type: ModelInfo
+    source = ""  # type: str
+    dtype = None  # type: np.dtype
+    fast = False  # type: bool
+    _program = None  # type: cl.Program
+    _kernels = None  # type: Dict[str, cl.Kernel]
     def __init__(self, source, model_info, dtype=generate.F32, fast=False):
         # type: (Dict[str,str], ModelInfo, np.dtype, bool) -> None
 …
         self.dtype = dtype
         self.fast = fast
-        self.timestamp = generate.ocl_timestamp(self.info)
-        self._cache_key = unique_id()
     def __getstate__(self):
 …
         # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None
         self.info, self.source, self.dtype, self.fast = state
+        self._program = self._kernels = None
     def make_kernel(self, q_vectors):
 …
         return GpuKernel(self, q_vectors)
+    @property
+    def Iq(self):
+        return self._fetch_kernel('Iq')
+    def fetch_kernel(self, name):
+    def get_function(self, name):
         # type: (str) -> cl.Kernel
         """
 …
         does not already exist.
         """
+        gpu = environment()
+        key = self._cache_key
+        if key not in gpu.cache:
+            program = gpu.compile_program(
+                self.info.name,
+                self.source['opencl'],
+                self.dtype,
+                self.fast,
+                self.timestamp)
+            variants = ['Iq', 'Iqxy', 'Imagnetic']
+            names = [generate.kernel_name(self.info, k) for k in variants]
+            kernels = [getattr(program, k) for k in names]
+            data = dict((k, v) for k, v in zip(variants, kernels))
+            # keep a handle to program so GC doesn't collect
+            data['program'] = program
+            gpu.cache[key] = data
+        else:
+            data = gpu.cache[key]
+        return data[name]
+# TODO: check that we don't need a destructor for buffers which go out of scope
+        if self._program is None:
+            self._prepare_program()
+        return self._kernels[name]
+    def _prepare_program(self):
+        # type: (str) -> None
+        env = environment()
+        timestamp = generate.ocl_timestamp(self.info)
+        program = env.compile_program(
+            self.info.name,
+            self.source['opencl'],
+            self.dtype,
+            self.fast,
+            timestamp)
+        variants = ['Iq', 'Iqxy', 'Imagnetic']
+        names = [generate.kernel_name(self.info, k) for k in variants]
+        functions = [getattr(program, k) for k in names]
+        self._kernels = {k: v for k, v in zip(variants, functions)}
+        # Keep a handle to program so GC doesn't collect.
+        self._program = program
+# TODO: Check that we don't need a destructor for buffers which go out of scope.
 class GpuInput(object):
     """
 …
     def __init__(self, q_vectors, dtype=generate.F32):
         # type: (List[np.ndarray], np.dtype) -> None
         # TODO: do we ever need double precision q?
+        # TODO: Do we ever need double precision q?
         self.nq = q_vectors[0].size
         self.dtype = np.dtype(dtype)
         self.is_2d = (len(q_vectors) == 2)
         # TODO: stretch input based on get_warp()
         # not doing it now since warp depends on kernel, which is not known
+        # TODO: Stretch input based on get_warp().
+        # Not doing it now since warp depends on kernel, which is not known
         # at this point, so instead using 32, which is good on the set of
         # architectures tested so far.
 …
             self.q[:self.nq] = q_vectors[0]
         self.global_size = [self.q.shape[0]]
+        self._cache_key = unique_id()
+    @property
+    def q_b(self):
+        """Lazy creation of q buffer so it can survive context reset"""
+        #print("creating inputs of size", self.global_size)
+        # Transfer input value to GPU.
         env = environment()
+        key = self._cache_key
+        if key not in env.cache:
+            context = env.context[self.dtype]
+            #print("creating inputs of size", self.global_size)
+            buffer = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                               hostbuf=self.q)
+            env.cache[key] = buffer
+        return env.cache[key]
+        context = env.context[self.dtype]
+        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                             hostbuf=self.q)
     def release(self):
         # type: () -> None
         """
+        Free the buffer associated with the q value
+        """
+        environment().free_buffer(id(self))
+        Free the buffer associated with the q value.
+        """
+        if self.q_b is not None:
+            self.q_b.release()
+            self.q_b = None
     def __del__(self):
 …
         self.release()
 class GpuKernel(Kernel):
     """
 …
     *model* is the GpuModel object to call
+    The following attributes are defined:
+    *info* is the module information
+    *dtype* is the kernel precision
+    *dim* is '1d' or '2d'
+    *result* is a vector to contain the results of the call
+    The resulting call method takes the *pars*, a list of values for
+    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
+    vectors for the polydisperse parameters.  *cutoff* determines the
+    integration limits: any points with combined weight less than *cutoff*
+    will not be calculated.
+    The kernel is derived from :class:`Kernel`, providing the
+    :meth:`call_kernel` method to evaluate the kernel for a given set of
+    parameters.  Because of the need to move the q values to the GPU before
+    evaluation, the kernel is instantiated for a particular set of q vectors,
+    and can be called many times without transfering q each time.
     Call :meth:`release` when done with the kernel instance.
     """
+    #: SAS model information structure.
+    info = None  # type: ModelInfo
+    #: Kernel precision.
+    dtype = None  # type: np.dtype
+    #: Kernel dimensions (1d or 2d).
+    dim = ""  # type: str
+    #: Calculation results, updated after each call to :meth:`_call_kernel`.
+    result = None  # type: np.ndarray
     def __init__(self, model, q_vectors):
         # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
+        # type: (GpuModel, List[np.ndarray]) -> None
         dtype = model.dtype
         self.q_input = GpuInput(q_vectors, dtype)
         self._model = model
+        # F16 isn't sufficient, so don't support it
+        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32
+        self._cache_key = unique_id()
+        # attributes accessed from the outside
+        # Attributes accessed from the outside.
         self.dim = '2d' if self.q_input.is_2d else '1d'
         self.info = model.info
+        self.dtype = model.dtype
+        # holding place for the returned value
+        self.dtype = dtype
+        # Converter to translate input to target type.
+        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32
+        # Holding place for the returned value.
         nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
+        extra_q = 4  # total weight, form volume, shell volume and R_eff
+        self.result = np.empty(self.q_input.nq*nout+extra_q, dtype)
+    @property
+    def _result_b(self):
+        """Lazy creation of result buffer so it can survive context reset"""
+        extra_q = 4  # Total weight, form volume, shell volume and R_eff.
+        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype)
+        # Allocate result value on GPU.
         env = environment()
+        key = self._cache_key
+        if key not in env.cache:
+            context = env.context[self.dtype]
+            width = ((self.result.size+31)//32)*32 * self.dtype.itemsize
+            buffer = cl.Buffer(context, mf.READ_WRITE, width)
+            env.cache[key] = buffer
+        return env.cache[key]
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        context = env.context[self.dtype]
+        width = ((self.result.size+31)//32)*32 * self.dtype.itemsize
+        self._result_b = cl.Buffer(context, mf.READ_WRITE, width)
+    def _call_kernel(self, call_details, values, cutoff, magnetic,
+                     radius_effective_mode):
+        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray
         env = environment()
         queue = env.queue[self._model.dtype]
         context = queue.context
+        # Arrange data transfer to/from card
+        q_b = self.q_input.q_b
+        result_b = self._result_b
+        # Arrange data transfer to card.
         details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
                               hostbuf=call_details.buffer)
 …
                              hostbuf=values)
+        # Setup kernel function and arguments.
         name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy'
         kernel = self._model.fetch_kernel(name)
+        kernel = self._model.get_function(name)
         kernel_args = [
+            np.uint32(self.q_input.nq), None, None,
+            details_b, values_b, q_b, result_b,
+            self._as_dtype(cutoff),
+            np.uint32(effective_radius_type),
+            np.uint32(self.q_input.nq),  # Number of inputs.
+            None,  # Placeholder for pd_start.
+            None,  # Placeholder for pd_stop.
+            details_b,  # Problem definition.
+            values_b,  # Parameter values.
+            self.q_input.q_b,  # Q values.
+            self._result_b,   # Result storage.
+            self._as_dtype(cutoff),  # Probability cutoff.
+            np.uint32(radius_effective_mode),  # R_eff mode.
+        ]
+        # Call kernel and retrieve results.
         #print("Calling OpenCL")
         #call_details.show(values)
-        #Call kernel and retrieve results
         wait_for = None
         last_nap = clock()
 …
                                *kernel_args, wait_for=wait_for)]
             if stop < call_details.num_eval:
                 # Allow other processes to run
+                # Allow other processes to run.
                 wait_for[0].wait()
                 current_time = clock()
 …
                     time.sleep(0.001)
                     last_nap = current_time
         cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for)
+        cl.enqueue_copy(queue, self.result, self._result_b, wait_for=wait_for)
         #print("result", self.result)
+        # Free buffers
+        for v in (details_b, values_b):
+            if v is not None:
+                v.release()
+        # Free buffers.
+        details_b.release()
+        values_b.release()
     def release(self):
 …
         Release resources associated with the kernel.
         """
-        environment().free_buffer(id(self))
         self.q_input.release()
+        if self._result_b is not None:
+            self._result_b.release()
+            self._result_b = None
     def __del__(self):

sasmodels/kernelcuda.py

-                      rf872fd1
+                      ra34b811
 import os
-import warnings
 import logging
 import time
 import re
+import atexit
 import numpy as np  # type: ignore
 # Attempt to setup cuda. This may fail if the pycuda package is not
+# Attempt to setup CUDA. This may fail if the pycuda package is not
 # installed or if it is installed but there are no devices available.
 try:
 …
 MAX_LOOPS = 2048
 def use_cuda():
+    env = os.environ.get("SAS_OPENCL", "").lower()
+    return HAVE_CUDA and (env == "" or env.startswith("cuda"))
+    # type: None -> bool
+    """Returns True if CUDA is the default compute engine."""
+    sas_opencl = os.environ.get("SAS_OPENCL", "CUDA").lower()
+    return HAVE_CUDA and sas_opencl.startswith("cuda")
 ENV = None
 …
     ENV = GpuEnvironment() if use_cuda() else None
 def environment():
     # type: () -> "GpuEnvironment"
 …
         if not HAVE_CUDA:
             raise RuntimeError("CUDA startup failed with ***"
                             + CUDA_ERROR + "***; using C compiler instead")
+                               + CUDA_ERROR + "***; using C compiler instead")
         reset_environment()
         if ENV is None:
 …
     return ENV
+# PyTest is not freeing ENV, so make sure it gets freed.
+atexit.register(lambda: ENV.release() if ENV is not None else None)
 def has_type(dtype):
     # type: (np.dtype) -> bool
 …
     Return true if device supports the requested precision.
     """
     # Assume the nvidia card supports 32-bit and 64-bit floats.
     # TODO: check if pycuda support F16
+    # Assume the NVIDIA card supports 32-bit and 64-bit floats.
+    # TODO: Check if pycuda support F16.
     return dtype in (generate.F32, generate.F64)
 FUNCTION_PATTERN = re.compile(r"""^
   (?P<space>\s*)                   # initial space
   (?P<qualifiers>^(?:\s*\b\w+\b\s*)+) # one or more qualifiers before function
   (?P<function>\s*\b\w+\b\s*[(])      # function name plus open parens
+  (?P<space>\s*)                       # Initial space.
+  (?P<qualifiers>^(?:\s*\b\w+\b\s*)+)  # One or more qualifiers before function.
+  (?P<function>\s*\b\w+\b\s*[(])       # Function name plus open parens.
   """, re.VERBOSE|re.MULTILINE)
 …
   """, re.VERBOSE|re.MULTILINE)
 def _add_device_tag(match):
     # type: (None) -> str
     # Note: should be re.Match, but that isn't a simple type
+    # Note: Should be re.Match, but that isn't a simple type.
     """
     replace qualifiers with __device__ qualifiers if needed
 …
         return "".join((space, "__device__ ", qualifiers, function))
 def mark_device_functions(source):
     # type: (str) -> str
 …
     return FUNCTION_PATTERN.sub(_add_device_tag, source)
 def show_device_functions(source):
     # type: (str) -> str
 …
     """
     for match in FUNCTION_PATTERN.finditer(source):
+        print(match.group('qualifiers').replace('\n',r'\n'), match.group('function'), '(')
+        print(match.group('qualifiers').replace('\n', r'\n'),
+              match.group('function'), '(')
     return source
 def compile_model(source, dtype, fast=False):
 …
     #options = ['--verbose', '-E']
     options = ['--use_fast_math'] if fast else None
     program = SourceModule(source, no_extern_c=True, options=options) # include_dirs=[...]
+    program = SourceModule(source, no_extern_c=True, options=options) #, include_dirs=[...])
     #print("done with "+program)
 …
 # for now, this returns one device in the context
 # TODO: create a context that contains all devices on all platforms
+# For now, this returns one device in the context.
+# TODO: Create a context that contains all devices on all platforms.
 class GpuEnvironment(object):
     """
     GPU context, with possibly many devices, and one queue per device.
+    GPU context for CUDA.
     """
     context = None # type: cuda.Context
     def __init__(self, devnum=None):
         # type: (int) -> None
-        # Byte boundary for data alignment
-        #self.data_boundary = max(d.min_data_type_align_size
-        #                         for d in self.context.devices)
-        self.compiled = {}
         env = os.environ.get("SAS_OPENCL", "").lower()
         if devnum is None and env.startswith("cuda:"):
             devnum = int(env[5:])
         # Set the global context to the particular device number if one is
         # given, otherwise use the default context.  Perhaps this will be set
 …
             self.context = make_default_context()
+        ## Byte boundary for data alignment.
+        #self.data_boundary = max(d.min_data_type_align_size
+        #                         for d in self.context.devices)
+        # Cache for compiled programs, and for items in context.
+        self.compiled = {}
     def release(self):
+        """Free the CUDA device associated with this context."""
         if self.context is not None:
             self.context.pop()
 …
         Compile the program for the device in the given context.
         """
+        # Note: PyOpenCL caches based on md5 hash of source, options and device
+        # so we don't really need to cache things for ourselves.  I'll do so
+        # anyway just to save some data munging time.
+        # Note: PyCuda (probably) caches but I'll do so as well just to
+        # save some data munging time.
         tag = generate.tag_source(source)
         key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else ""))
         # Check timestamp on program
+        # Check timestamp on program.
         program, program_timestamp = self.compiled.get(key, (None, np.inf))
         if program_timestamp < timestamp:
 …
         return program
 class GpuModel(KernelModel):
     """
 …
     that the compiler is allowed to take shortcuts.
     """
     info = None # type: ModelInfo
     source = "" # type: str
     dtype = None # type: np.dtype
     fast = False # type: bool
     program = None # type: SourceModule
     _kernels = None # type: List[cuda.Function]
+    info = None  # type: ModelInfo
+    source = ""  # type: str
+    dtype = None  # type: np.dtype
+    fast = False  # type: bool
+    _program = None  # type: SourceModule
+    _kernels = None  # type: Dict[str, cuda.Function]
     def __init__(self, source, model_info, dtype=generate.F32, fast=False):
 …
         self.dtype = dtype
         self.fast = fast
-        self.program = None # delay program creation
-        self._kernels = None
     def __getstate__(self):
 …
         # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None
         self.info, self.source, self.dtype, self.fast = state
         self.program = None
+        self._program = self._kernels = None
     def make_kernel(self, q_vectors):
         # type: (List[np.ndarray]) -> "GpuKernel"
+        if self.program is None:
+            compile_program = environment().compile_program
+            timestamp = generate.ocl_timestamp(self.info)
+            self.program = compile_program(
+                self.info.name,
+                self.source['opencl'],
+                self.dtype,
+                self.fast,
+                timestamp)
+            variants = ['Iq', 'Iqxy', 'Imagnetic']
+            names = [generate.kernel_name(self.info, k) for k in variants]
+            kernels = [self.program.get_function(k) for k in names]
+            self._kernels = dict((k, v) for k, v in zip(variants, kernels))
+        is_2d = len(q_vectors) == 2
+        if is_2d:
+            kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']]
+        else:
+            kernel = [self._kernels['Iq']]*2
+        return GpuKernel(kernel, self.dtype, self.info, q_vectors)
+    def release(self):
+        # type: () -> None
+        """
+        Free the resources associated with the model.
+        """
+        if self.program is not None:
+            self.program = None
+    def __del__(self):
+        # type: () -> None
+        self.release()
+# TODO: check that we don't need a destructor for buffers which go out of scope
+        return GpuKernel(self, q_vectors)
+    def get_function(self, name):
+        # type: (str) -> cuda.Function
+        """
+        Fetch the kernel from the environment by name, compiling it if it
+        does not already exist.
+        """
+        if self._program is None:
+            self._prepare_program()
+        return self._kernels[name]
+    def _prepare_program(self):
+        # type: (str) -> None
+        env = environment()
+        timestamp = generate.ocl_timestamp(self.info)
+        program = env.compile_program(
+            self.info.name,
+            self.source['opencl'],
+            self.dtype,
+            self.fast,
+            timestamp)
+        variants = ['Iq', 'Iqxy', 'Imagnetic']
+        names = [generate.kernel_name(self.info, k) for k in variants]
+        functions = [program.get_function(k) for k in names]
+        self._kernels = {k: v for k, v in zip(variants, functions)}
+        # Keep a handle to program so GC doesn't collect.
+        self._program = program
+# TODO: Check that we don't need a destructor for buffers which go out of scope.
 class GpuInput(object):
     """
 …
     def __init__(self, q_vectors, dtype=generate.F32):
         # type: (List[np.ndarray], np.dtype) -> None
         # TODO: do we ever need double precision q?
+        # TODO: Do we ever need double precision q?
         self.nq = q_vectors[0].size
         self.dtype = np.dtype(dtype)
         self.is_2d = (len(q_vectors) == 2)
         # TODO: stretch input based on get_warp()
         # not doing it now since warp depends on kernel, which is not known
+        # TODO: Stretch input based on get_warp().
+        # Not doing it now since warp depends on kernel, which is not known
         # at this point, so instead using 32, which is good on the set of
         # architectures tested so far.
         if self.is_2d:
+            # Note: 16 rather than 15 because result is 1 longer than input.
+            width = ((self.nq+16)//16)*16
+            width = ((self.nq+15)//16)*16
             self.q = np.empty((width, 2), dtype=dtype)
             self.q[:self.nq, 0] = q_vectors[0]
             self.q[:self.nq, 1] = q_vectors[1]
         else:
+            # Note: 32 rather than 31 because result is 1 longer than input.
+            width = ((self.nq+32)//32)*32
+            width = ((self.nq+31)//32)*32
             self.q = np.empty(width, dtype=dtype)
             self.q[:self.nq] = q_vectors[0]
         self.global_size = [self.q.shape[0]]
         #print("creating inputs of size", self.global_size)
+        # Transfer input value to GPU.
         self.q_b = cuda.to_device(self.q)
 …
         # type: () -> None
         """
         Free the memory.
+        Free the buffer associated with the q value.
         """
         if self.q_b is not None:
 …
         self.release()
 class GpuKernel(Kernel):
     """
     Callable SAS kernel.
+    *kernel* is the GpuKernel object to call
+    *model_info* is the module information
+    *q_vectors* is the q vectors at which the kernel should be evaluated
+    *dtype* is the kernel precision
+    The resulting call method takes the *pars*, a list of values for
+    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
+    vectors for the polydisperse parameters.  *cutoff* determines the
+    integration limits: any points with combined weight less than *cutoff*
+    will not be calculated.
+    *model* is the GpuModel object to call
+    The kernel is derived from :class:`Kernel`, providing the
+    :meth:`call_kernel` method to evaluate the kernel for a given set of
+    parameters.  Because of the need to move the q values to the GPU before
+    evaluation, the kernel is instantiated for a particular set of q vectors,
+    and can be called many times without transfering q each time.
     Call :meth:`release` when done with the kernel instance.
     """
+    def __init__(self, kernel, dtype, model_info, q_vectors):
+        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
+    #: SAS model information structure.
+    info = None  # type: ModelInfo
+    #: Kernel precision.
+    dtype = None  # type: np.dtype
+    #: Kernel dimensions (1d or 2d).
+    dim = ""  # type: str
+    #: Calculation results, updated after each call to :meth:`_call_kernel`.
+    result = None  # type: np.ndarray
+    def __init__(self, model, q_vectors):
+        # type: (GpuModel, List[np.ndarray]) -> None
+        dtype = model.dtype
         self.q_input = GpuInput(q_vectors, dtype)
+        self.kernel = kernel
+        # F16 isn't sufficient, so don't support it
+        self._model = model
+        # Attributes accessed from the outside.
+        self.dim = '2d' if self.q_input.is_2d else '1d'
+        self.info = model.info
+        self.dtype = dtype
+        # Converter to translate input to target type.
         self._as_dtype = np.float64 if dtype == generate.F64 else np.float32
+        # attributes accessed from the outside
+        self.dim = '2d' if self.q_input.is_2d else '1d'
+        self.info = model_info
+        self.dtype = dtype
+        # holding place for the returned value
+        # Holding place for the returned value.
         nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
+        extra_q = 4  # total weight, form volume, shell volume and R_eff
+        self.result = np.empty(self.q_input.nq*nout+extra_q, dtype)
+        # Inputs and outputs for each kernel call
+        # Note: res may be shorter than res_b if global_size != nq
+        extra_q = 4  # Total weight, form volume, shell volume and R_eff.
+        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype)
+        # Allocate result value on GPU.
         width = ((self.result.size+31)//32)*32 * self.dtype.itemsize
+        self.result_b = cuda.mem_alloc(width)
+        self._need_release = [self.result_b]
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
+        # Arrange data transfer to card
+        self._result_b = cuda.mem_alloc(width)
+    def _call_kernel(self, call_details, values, cutoff, magnetic,
+                     radius_effective_mode):
+        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray
+        # Arrange data transfer to card.
         details_b = cuda.to_device(call_details.buffer)
         values_b = cuda.to_device(values)
+        kernel = self.kernel[1 if magnetic else 0]
+        args = [
+            np.uint32(self.q_input.nq), None, None,
+            details_b, values_b, self.q_input.q_b, self.result_b,
+            self._as_dtype(cutoff),
+            np.uint32(effective_radius_type),
+        # Setup kernel function and arguments.
+        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy'
+        kernel = self._model.get_function(name)
+        kernel_args = [
+            np.uint32(self.q_input.nq),  # Number of inputs.
+            None,  # Placeholder for pd_start.
+            None,  # Placeholder for pd_stop.
+            details_b,  # Problem definition.
+            values_b,  # Parameter values.
+            self.q_input.q_b,  # Q values.
+            self._result_b,   # Result storage.
+            self._as_dtype(cutoff),  # Probability cutoff.
+            np.uint32(radius_effective_mode),  # R_eff mode.
+        ]
         grid = partition(self.q_input.nq)
+        #print("Calling OpenCL")
+        # Call kernel and retrieve results.
+        #print("Calling CUDA")
         #call_details.show(values)
-        # Call kernel and retrieve results
         last_nap = time.clock()
         step = 100000000//self.q_input.nq + 1
 …
             stop = min(start + step, call_details.num_eval)
             #print("queuing",start,stop)
             args[1:3] = [np.int32(start), np.int32(stop)]
             kernel(*args, **grid)
+            kernel_args[1:3] = [np.int32(start), np.int32(stop)]
+            kernel(*kernel_args, **grid)
             if stop < call_details.num_eval:
                 sync()
                 # Allow other processes to run
+                # Allow other processes to run.
                 current_time = time.clock()
                 if current_time - last_nap > 0.5:
 …
                     last_nap = current_time
         sync()
         cuda.memcpy_dtoh(self.result, self.result_b)
+        cuda.memcpy_dtoh(self.result, self._result_b)
         #print("result", self.result)
 …
         Release resources associated with the kernel.
         """
+        for p in self._need_release:
+            p.free()
+        self._need_release = []
+        self.q_input.release()
+        if self._result_b is not None:
+            self._result_b.free()
+            self._result_b = None
     def __del__(self):
 …
     Note: Maybe context.synchronize() is sufficient.
     """
+    #return # The following works in C++; don't know what pycuda is doing
+    # Create an event with which to synchronize
+    # Create an event with which to synchronize.
     done = cuda.Event()
 …
     done.record()
     #line added to not hog resources
+    # Make sure we don't hog resource while waiting to sync.
     while not done.query():
         time.sleep(0.01)
 …
     # Block until the GPU executes the kernel.
     done.synchronize()
     # Clean up the event; I don't think they can be reused.
     del done

sasmodels/kerneldll.py

-                      re44432d
+                      ra34b811
 # pylint: enable=unused-import
 # Compiler output is a byte stream that needs to be decode in python 3
+# Compiler output is a byte stream that needs to be decode in python 3.
 decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8'))
 …
         COMPILER = "tinycc"
     elif "VCINSTALLDIR" in os.environ:
+        # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment
+        # and we can use the MSVC compiler.  Otherwise, if tinycc is available
+        # the use it.  Otherwise, hope that mingw is available.
+        # If vcvarsall.bat has been called, then VCINSTALLDIR is in the
+        # environment and we can use the MSVC compiler.  Otherwise, if
+        # tinycc is available then use it.  Otherwise, hope that mingw
+        # is available.
         COMPILER = "msvc"
     else:
 …
     COMPILER = "unix"
 ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86
+ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86.
 if COMPILER == "unix":
     # Generic unix compile
     # On mac users will need the X code command line tools installed
+    # Generic unix compile.
+    # On Mac users will need the X code command line tools installed.
     #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
     CC = "cc -shared -fPIC -std=c99 -O2 -Wall".split()
     # add openmp support if not running on a mac
+    # Add OpenMP support if not running on a Mac.
     if sys.platform != "darwin":
         # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9)
+        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9).
         # Shut it off for all unix until we can investigate.
         #CC.append("-fopenmp")
 …
     # vcomp90.dll on the path.  One may be found here:
     #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll
     # Copy this to the python directory and uncomment the OpenMP COMPILE
     # TODO: remove intermediate OBJ file created in the directory
     # TODO: maybe don't use randomized name for the c file
     # TODO: maybe ask distutils to find MSVC
+    # Copy this to the python directory and uncomment the OpenMP COMPILE.
+    # TODO: Remove intermediate OBJ file created in the directory.
+    # TODO: Maybe don't use randomized name for the c file.
+    # TODO: Maybe ask distutils to find MSVC.
     CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split()
     if "SAS_OPENMP" in os.environ:
 …
 ALLOW_SINGLE_PRECISION_DLLS = True
+def compile(source, output):
+def compile_model(source, output):
     # type: (str, str) -> None
     """
 …
     logging.info(command_str)
     try:
         # need shell=True on windows to keep console box from popping up
+        # Need shell=True on windows to keep console box from popping up.
         shell = (os.name == 'nt')
         subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT)
 …
         raise RuntimeError("compile failed.  File is in %r"%source)
 def dll_name(model_info, dtype):
     # type: (ModelInfo, np.dtype) ->  str
 …
     basename += ARCH + ".so"
     # Hack to find precompiled dlls
+    # Hack to find precompiled dlls.
     path = joinpath(generate.DATA_PATH, '..', 'compiled_models', basename)
     if os.path.exists(path):
 …
         raise ValueError("16 bit floats not supported")
     if dtype == F32 and not ALLOW_SINGLE_PRECISION_DLLS:
         dtype = F64  # Force 64-bit dll
     # Note: dtype may be F128 for long double precision
+        dtype = F64  # Force 64-bit dll.
+    # Note: dtype may be F128 for long double precision.
     dll = dll_path(model_info, dtype)
 …
         need_recompile = dll_time < newest_source
     if need_recompile:
         # Make sure the DLL path exists
+        # Make sure the DLL path exists.
         if not os.path.exists(SAS_DLL_PATH):
             os.makedirs(SAS_DLL_PATH)
 …
         with os.fdopen(system_fd, "w") as file_handle:
             file_handle.write(source)
         compile(source=filename, output=dll)
         # comment the following to keep the generated c file
         # Note: if there is a syntax error then compile raises an error
+        compile_model(source=filename, output=dll)
+        # Comment the following to keep the generated C file.
+        # Note: If there is a syntax error then compile raises an error
         # and the source file will not be deleted.
         os.unlink(filename)
 …
         self.dllpath = dllpath
         self._dll = None  # type: ct.CDLL
         self._kernels = None # type: List[Callable, Callable]
+        self._kernels = None  # type: List[Callable, Callable]
         self.dtype = np.dtype(dtype)
 …
         # type: (List[np.ndarray]) -> DllKernel
         q_input = PyInput(q_vectors, self.dtype)
         # Note: pickle not supported for DllKernel
+        # Note: DLL is lazy loaded.
         if self._dll is None:
             self._load_dll()
 …
         self._dll = None
 class DllKernel(Kernel):
     """
 …
     def __init__(self, kernel, model_info, q_input):
         # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None
+        #,model_info,q_input)
+        dtype = q_input.dtype
+        self.q_input = q_input
         self.kernel = kernel
+        # Attributes accessed from the outside.
+        self.dim = '2d' if q_input.is_2d else '1d'
         self.info = model_info
+        self.q_input = q_input
+        self.dtype = q_input.dtype
+        self.dim = '2d' if q_input.is_2d else '1d'
+        # leave room for f1/f2 results in case we need to compute beta for 1d models
+        self.dtype = dtype
+        # Converter to translate input to target type.
+        self._as_dtype = (np.float32 if dtype == generate.F32
+                          else np.float64 if dtype == generate.F64
+                          else np.float128)
+        # Holding place for the returned value.
         nout = 2 if self.info.have_Fq else 1
         # +4 for total weight, shell volume, effective radius, form volume
         self.result = np.empty(q_input.nq*nout + 4, self.dtype)
+        self.real = (np.float32 if self.q_input.dtype == generate.F32
                      else np.float64 if self.q_input.dtype == generate.F64
                      else np.float128)
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        extra_q = 4  # Total weight, form volume, shell volume and R_eff.
+        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype)
+    def _call_kernel(self, call_details, values, cutoff, magnetic,
+                     radius_effective_mode):
+        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray
+        # Setup kernel function and arguments.
         kernel = self.kernel[1 if magnetic else 0]
         args = [
             self.q_input.nq, # nq
             None, # pd_start
             None, # pd_stop pd_stride[MAX_PD]
             call_details.buffer.ctypes.data, # problem
             values.ctypes.data,  # pars
             self.q_input.q.ctypes.data, # q
             self.result.ctypes.data,   # results
             self.real(cutoff), # cutoff
             effective_radius_type, # cutoff
+        kernel_args = [
+            self.q_input.nq,  # Number of inputs.
+            None,  # Placeholder for pd_start.
+            None,  # Placeholder for pd_stop.
+            call_details.buffer.ctypes.data,  # Problem definition.
+            values.ctypes.data,  # Parameter values.
+            self.q_input.q.ctypes.data,  # Q values.
+            self.result.ctypes.data,   # Result storage.
+            self._as_dtype(cutoff),  # Probability cutoff.
+            radius_effective_mode,  # R_eff mode.
+        ]
+        # Call kernel and retrieve results.
         #print("Calling DLL")
         #call_details.show(values)
         step = 100
+        # TODO: Do we need the explicit sleep like the OpenCL and CUDA loops?
         for start in range(0, call_details.num_eval, step):
             stop = min(start + step, call_details.num_eval)
             args[1:3] = [start, stop]
             kernel(*args) # type: ignore
+            kernel_args[1:3] = [start, stop]
+            kernel(*kernel_args) # type: ignore
     def release(self):
         # type: () -> None
         """
         Release any resources associated with the kernel.
+        Release resources associated with the kernel.
         """
+        self.q_input.release()
+        # TODO: OpenCL/CUDA allocate q_input in __init__ and free it in release.
+        # Should we be doing the same for DLL?
+        #self.q_input.release()
+        pass
+    def __del__(self):
+        # type: () -> None
+        self.release()

sasmodels/kernelpy.py

-                      raa8c6e0
+                      ra34b811
     from numpy import cbrt
 except ImportError:
+    def cbrt(x): return x ** (1.0/3.0)
+    def cbrt(x):
+        """Return cubed root of x."""
+        return x ** (1.0/3.0)
 from .generate import F64
 …
 logger = logging.getLogger(__name__)
 class PyModel(KernelModel):
     """
 …
     """
     def __init__(self, model_info):
         # Make sure Iq is available and vectorized
+        # Make sure Iq is available and vectorized.
         _create_default_functions(model_info)
         self.info = model_info
         self.dtype = np.dtype('d')
         logger.info("make python model " + self.info.name)
+        logger.info("make python model %s", self.info.name)
     def make_kernel(self, q_vectors):
+        """Instantiate the python kernel with input *q_vectors*"""
         q_input = PyInput(q_vectors, dtype=F64)
         return PyKernel(self.info, q_input)
 …
         """
         pass
 class PyInput(object):
 …
         self.q = None
 class PyKernel(Kernel):
     """
 …
         parameter_vector = np.empty(len(partable.call_parameters)-2, 'd')
         # Create views into the array to hold the arguments
+        # Create views into the array to hold the arguments.
         offset = 0
         kernel_args, volume_args = [], []
 …
         volume = model_info.form_volume
         shell = model_info.shell_volume
         radius = model_info.effective_radius
+        radius = model_info.radius_effective
         self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume
                         else (lambda: [volume(*volume_args)]*2) if volume
 …
                         else (lambda mode: 1.0))
+    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
+    def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_mode):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
         if magnetic:
 …
         #print("Calling python kernel")
         #call_details.show(values)
         radius = ((lambda: 0.0) if effective_radius_type == 0
                   else (lambda: self._radius(effective_radius_type)))
+        radius = ((lambda: 0.0) if radius_effective_mode == 0
+                  else (lambda: self._radius(radius_effective_mode)))
         self.result = _loops(
             self._parameter_vector, self._form, self._volume, radius,
 …
         self.q_input.release()
         self.q_input = None
 def _loops(parameters,    # type: np.ndarray
 …
         total = np.zeros(nq, 'd')
         for loop_index in range(call_details.num_eval):
             # update polydispersity parameter values
+            # Update polydispersity parameter values.
             if p0_index == p0_length:
                 pd_index = (loop_index//pd_stride)%pd_length
 …
             p0_index += 1
             if weight > cutoff:
                 # Call the scattering function
+                # Call the scattering function.
                 # Assume that NaNs are only generated if the parameters are bad;
                 # exclude all q for that NaN.  Even better would be to have an
 …
                     continue
                 # update value and norm
+                # Update value and norm.
                 total += weight * Iq
                 weight_norm += weight
 …
     any functions that are not already marked as vectorized.
     """
     # Note: must call create_vector_Iq before create_vector_Iqxy
+    # Note: Must call create_vector_Iq before create_vector_Iqxy.
     _create_vector_Iq(model_info)
     _create_vector_Iqxy(model_info)

sasmodels/list_pars.py

-                      r2d81cfe
+                      rb297ba9
 from .compare import columnize
 def find_pars(type=None):
+def find_pars(kind=None):
     """
     Find all parameters in all models.
 …
         model_info = load_model_info(name)
         for p in model_info.parameters.kernel_parameters:
             if type is None or p.type == type:
+            if kind is None or p.type == kind:
                 partable.setdefault(p.name, [])
                 partable[p.name].append(name)
     return partable
 def list_pars(names_only=True, type=None):
+def list_pars(names_only=True, kind=None):
     """
     Print all parameters in all models.
 …
     occurs in.
     """
     partable = find_pars(type)
+    partable = find_pars(kind)
     if names_only:
         print(columnize(list(sorted(partable.keys()))))
 …
         help="list models which use this argument")
     parser.add_argument(
         'type', default="any", nargs='?',
+        'kind', 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")
+        type=lambda v: None if v == 'any' else None if v == 'none' else v,
+        help="only list arguments of the given kind")
     args = parser.parse_args()
     list_pars(names_only=not args.verbose, type=args.type)
+    list_pars(names_only=not args.verbose, kind=args.kind)
 if __name__ == "__main__":

sasmodels/mixture.py

-                      r39a06c9
+                      rb297ba9
     def random():
+        """Random set of model parameters for mixture model"""
         combined_pars = {}
         for k, part in enumerate(parts):
 …
 class MixtureModel(KernelModel):
+    """
+    Model definition for mixture of models.
+    """
     def __init__(self, model_info, parts):
         # type: (ModelInfo, List[KernelModel]) -> None
 …
         kernels = [part.make_kernel(q_vectors) for part in self.parts]
         return MixtureKernel(self.info, kernels)
+    make_kernel.__doc__ = KernelModel.make_kernel.__doc__
     def release(self):
         # type: () -> None
+        """
+        Free resources associated with the model.
+        """
+        """Free resources associated with the model."""
         for part in self.parts:
             part.release()
+    release.__doc__ = KernelModel.release.__doc__
 class MixtureKernel(Kernel):
+    """
+    Instantiated kernel for mixture of models.
+    """
     def __init__(self, model_info, kernels):
         # type: (ModelInfo, List[Kernel]) -> None
 …
         self.results = []  # type: List[np.ndarray]
     def __call__(self, call_details, values, cutoff, magnetic):
+    def Iq(self, call_details, values, cutoff, magnetic):
         # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray
         scale, background = values[0:2]
 …
         # remember the parts for plotting later
         self.results = []  # type: List[np.ndarray]
         parts = MixtureParts(self.info, self.kernels, call_details, values)
+        parts = _MixtureParts(self.info, self.kernels, call_details, values)
         for kernel, kernel_details, kernel_values in parts:
             #print("calling kernel", kernel.info.name)
 …
         return scale*total + background
+    Iq.__doc__ = Kernel.Iq.__doc__
+    __call__ = Iq
     def release(self):
         # type: () -> None
+        """Free resources associated with the kernel."""
         for k in self.kernels:
             k.release()
+class MixtureParts(object):
+# Note: _MixtureParts doesn't implement iteration correctly, and only allows
+# a single iterator to be active at once.  It doesn't matter in this case
+# since _MixtureParts is only used in one place, but it is not clean style.
+class _MixtureParts(object):
+    """
+    Mixture component iterator.
+    """
     def __init__(self, model_info, kernels, call_details, values):
         # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None
 …
         self.values = values
         self.spin_index = model_info.parameters.npars + 2
+        # The following are redefined by __iter__, but set them here so that
+        # lint complains a little less.
+        self.part_num = -1
+        self.par_index = -1
+        self.mag_index = -1
         #call_details.show(values)

sasmodels/model_test.py

-                      r5024a56
+                      r8795b6f
 Usage::
+    python -m sasmodels.model_test [opencl|cuda|dll] model1 model2 ...
+    if model1 is 'all', then all except the remaining models will be tested
+    python -m sasmodels.model_test [opencl|cuda|dll|all] model1 model2 ...
+If model1 is 'all', then all except the remaining models will be tested.
+Subgroups are also possible, such as 'py', 'single' or '1d'.  See
+:func:`core.list_models` for details.
 Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1),
 …
 from __future__ import print_function
+import argparse
 import sys
 import unittest
 …
 import numpy as np  # type: ignore
-from . import core
 from .core import list_models, load_model_info, build_model
 from .direct_model import call_kernel, call_Fq
 …
     suite = unittest.TestSuite()
+    if models[0] in core.KINDS:
+    try:
+        # See if the first model parses as a model group
+        group = list_models(models[0])
         skip = models[1:]
         models = list_models(models[0])
     else:
+        models = group
+    except Exception:
         skip = []
     for model_name in models:
 …
         test_method_name = "test_%s_python" % model_info.id
         test = ModelTestCase(test_name, model_info,
                                 test_method_name,
                                 platform="dll",  # so that
                                 dtype="double",
                                 stash=stash)
+                             test_method_name,
+                             platform="dll",  # so that
+                             dtype="double",
+                             stash=stash)
         suite.addTest(test)
     else:   # kernel implemented in C
 …
             test_method_name = "test_%s_dll" % model_info.id
             test = ModelTestCase(test_name, model_info,
                                     test_method_name,
                                     platform="dll",
                                     dtype="double",
                                     stash=stash)
+                                 test_method_name,
+                                 platform="dll",
+                                 dtype="double",
+                                 stash=stash)
             suite.addTest(test)
 …
             # presence of *single=False* in the model file.
             test = ModelTestCase(test_name, model_info,
                                     test_method_name,
                                     platform="ocl", dtype=None,
                                     stash=stash)
+                                 test_method_name,
+                                 platform="ocl", dtype=None,
+                                 stash=stash)
             #print("defining", test_name)
             suite.addTest(test)
 …
         # test using cuda if desired and available
         if 'cuda' in loaders and use_cuda():
             test_name = "%s-cuda"%model_name
+            test_name = "%s-cuda" % model_info.id
             test_method_name = "test_%s_cuda" % model_info.id
             # Using dtype=None so that the models that are only
 …
             # presence of *single=False* in the model file.
             test = ModelTestCase(test_name, model_info,
                                     test_method_name,
                                     platform="cuda", dtype=None,
                                     stash=stash)
+                                 test_method_name,
+                                 platform="cuda", dtype=None,
+                                 stash=stash)
             #print("defining", test_name)
             suite.addTest(test)
 …
                     s_name = pars.pop('@S')
                     ps_test = [pars] + list(test[1:])
+                    #print("PS TEST PARAMS!!!",ps_test)
                     # build the P@S model
                     s_info = load_model_info(s_name)
 …
                                            platform=self.platform)
                     # run the tests
+                    #self.info = ps_model.info
+                    #print("SELF.INFO PARAMS!!!",[p.id for p in self.info.parameters.call_parameters])
+                    #print("PS MODEL PARAMETERS:",[p.id for p in ps_model.info.parameters.call_parameters])
                     results.append(self.run_one(ps_model, ps_test))
 …
             """Run a single test case."""
             user_pars, x, y = test[:3]
+            pars = expand_pars(self.info.parameters, user_pars)
+            invalid = invalid_pars(self.info.parameters, pars)
+            #print("PS MODEL PARAMETERS:",[p.id for p in model.info.parameters.call_parameters])
+            pars = expand_pars(model.info.parameters, user_pars)
+            invalid = invalid_pars(model.info.parameters, pars)
             if invalid:
                 raise ValueError("Unknown parameters in test: " + ", ".join(invalid))
 …
             else:
                 y1 = y
                 y2 = test[3] if not isinstance(test[3], list) else [test[3]]
                 F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars)
                 if F1 is not None:  # F1 is none for models with Iq instead of Fq
                     self._check_vectors(x, y1, F1, 'F')
                 self._check_vectors(x, y2, F2, 'F^2')
+                y2 = test[3] if isinstance(test[3], list) else [test[3]]
+                F, Fsq, R_eff, volume, volume_ratio = call_Fq(kernel, pars)
+                if F is not None:  # F is none for models with Iq instead of Fq
+                    self._check_vectors(x, y1, F, 'F')
+                self._check_vectors(x, y2, Fsq, 'F^2')
                 self._check_scalar(test[4], R_eff, 'R_eff')
                 self._check_scalar(test[5], volume, 'volume')
                 self._check_scalar(test[6], volume_ratio, 'form:shell ratio')
                 return F2
+                return Fsq
         def _check_scalar(self, target, actual, name):
+            if target is None:
+                # smoke test --- make sure it runs and produces a value
+                self.assertTrue(not np.isnan(actual),
+                                'invalid %s: %s' % (name, actual))
+            elif np.isnan(target):
+                # make sure nans match
+                self.assertTrue(np.isnan(actual),
+                                '%s: expected:%s; actual:%s'
+                                % (name, target, actual))
+            else:
+                # is_near does not work for infinite values, so also test
+                # for exact values.
+                self.assertTrue(target == actual or is_near(target, actual, 5),
+                                '%s: expected:%s; actual:%s'
+                                % (name, target, actual))
+            self.assertTrue(is_near(target, actual, 5),
+                            '%s: expected:%s; actual:%s'
+                            % (name, target, actual))
         def _check_vectors(self, x, target, actual, name='I'):
 …
                              '%s(...) returned wrong length'%name)
             for xi, yi, actual_yi in zip(x, target, actual):
+                if yi is None:
+                    # smoke test --- make sure it runs and produces a value
+                    self.assertTrue(not np.isnan(actual_yi),
+                                    'invalid %s(%s): %s' % (name, xi, actual_yi))
+                elif np.isnan(yi):
+                    # make sure nans match
+                    self.assertTrue(np.isnan(actual_yi),
+                                    '%s(%s): expected:%s; actual:%s'
+                                    % (name, xi, yi, actual_yi))
+                else:
+                    # is_near does not work for infinite values, so also test
+                    # for exact values.
+                    self.assertTrue(yi == actual_yi or is_near(yi, actual_yi, 5),
+                                    '%s(%s); expected:%s; actual:%s'
+                                    % (name, xi, yi, actual_yi))
+                self.assertTrue(is_near(yi, actual_yi, 5),
+                                '%s(%s): expected:%s; actual:%s'
+                                % (name, xi, target, actual))
     return ModelTestCase
 …
     invalid = []
     for par in sorted(pars.keys()):
+        # special handling of R_eff mode, which is not a usual parameter
+        # Ignore the R_eff mode parameter when checking for valid parameters.
+        # It is an allowed parameter for a model even though it does not exist
+        # in the parameter table.  The call_Fq() function pops it from the
+        # parameter list and sends it directly to kernel.Fq().
         if par == product.RADIUS_MODE_ID:
             continue
 …
     """
     Returns true if *actual* is within *digits* significant digits of *target*.
+    """
+    import math
+    shift = 10**math.ceil(math.log10(abs(target)))
+    return abs(target-actual)/shift < 1.5*10**-digits
+    *taget* zero and inf should match *actual* zero and inf.  If you want to
+    accept eps for zero, choose a value such as 1e-10, which must match up to
+    +/- 1e-15 when *digits* is the default value of 5.
+    If *target* is None, then just make sure that *actual* is not NaN.
+    If *target* is NaN, make sure *actual* is NaN.
+    """
+    if target is None:
+        # target is None => actual cannot be NaN
+        return not np.isnan(actual)
+    elif target == 0.:
+        # target is 0. => actual must be 0.
+        # Note: if small values are allowed, then use maybe test zero against eps instead?
+        return actual == 0.
+    elif np.isfinite(target):
+        shift = np.ceil(np.log10(abs(target)))
+        return abs(target-actual) < 1.5*10**(shift-digits)
+    elif target == actual:
+        # target is inf => actual must be inf of same sign
+        return True
+    else:
+        # target is NaN => actual must be NaN
+        return np.isnan(target) == np.isnan(actual)
 # CRUFT: old interface; should be deprecated and removed
 def run_one(model_name):
+    # type: (str) -> str
+    """
+    [Deprecated] Run the tests associated with *model_name*.
+    Use the following instead::
+        succss, output = check_model(load_model_info(model_name))
+    """
     # msg = "use check_model(model_info) rather than run_one(model_name)"
     # warnings.warn(msg, category=DeprecationWarning, stacklevel=2)
 …
         return output
     success, output = check_model(model_info)
+    _, output = check_model(model_info)
     return output
 def check_model(model_info):
     # type: (ModelInfo) -> str
+    # type: (ModelInfo) -> Tuple[bool, str]
     """
     Run the tests for a single model, capturing the output.
 …
-def main(*models):
-    # type: (*str) -> int
-    """
-    Run tests given is models.
-    Returns 0 if success or 1 if any tests fail.
-    """
-    try:
-        from xmlrunner import XMLTestRunner as TestRunner
-        test_args = {'output': 'logs'}
-    except ImportError:
-        from unittest import TextTestRunner as TestRunner
-        test_args = {}
-    if models and models[0] == '-v':
-        verbosity = 2
-        models = models[1:]
-    else:
-        verbosity = 1
-    if models and models[0] == 'opencl':
-        if not use_opencl():
-            print("opencl is not available")
-            return 1
-        loaders = ['opencl']
-        models = models[1:]
-    elif models and models[0] == 'cuda':
-        if not use_cuda():
-            print("cuda is not available")
-            return 1
-        loaders = ['cuda']
-        models = models[1:]
-    elif models and models[0] == 'dll':
-        # TODO: test if compiler is available?
-        loaders = ['dll']
-        models = models[1:]
-    else:
-        loaders = ['dll']
-        if use_opencl():
-            loaders.append('opencl')
-        if use_cuda():
-            loaders.append('cuda')
-    if not models:
-        print("""\
-usage:
-  python -m sasmodels.model_test [-v] [opencl|cuda|dll] model1 model2 ...
-If -v is included on the command line, then use verbose output.
-If no platform is specified, then models will be tested with dll, and
-if available, OpenCL and CUDA; the compute target is ignored for pure python models.
-If model1 is 'all', then all except the remaining models will be tested.
-""")
-        return 1
-    runner = TestRunner(verbosity=verbosity, **test_args)
-    result = runner.run(make_suite(loaders, models))
-    return 1 if result.failures or result.errors else 0
 def model_tests():
     # type: () -> Iterator[Callable[[], None]]
 …
         loaders.append('cuda')
     tests = make_suite(loaders, ['all'])
     def build_test(test):
+    def _build_test(test):
         # In order for nosetest to show the test name, wrap the test.run_all
         # instance in function that takes the test name as a parameter which
 …
     for test in tests:
+        yield build_test(test)
+        yield _build_test(test)
+def main():
+    # type: (*str) -> int
+    """
+    Run tests given is models.
+    Returns 0 if success or 1 if any tests fail.
+    """
+    try:
+        from xmlrunner import XMLTestRunner as TestRunner
+        test_args = {'output': 'logs'}
+    except ImportError:
+        from unittest import TextTestRunner as TestRunner
+        test_args = {}
+    parser = argparse.ArgumentParser(description="Test SasModels Models")
+    parser.add_argument("-v", "--verbose", action="store_const",
+                        default=1, const=2, help="Use verbose output")
+    parser.add_argument("-e", "--engine", default="all",
+                        help="Engines on which to run the test.  "
+                        "Valid values are opencl, cuda, dll, and all. "
+                        "Defaults to all if no value is given")
+    parser.add_argument("models", nargs="*",
+                        help="The names of the models to be tested.  "
+                        "If the first model is 'all', then all but the listed "
+                        "models will be tested.  See core.list_models() for "
+                        "names of other groups, such as 'py' or 'single'.")
+    opts = parser.parse_args()
+    if opts.engine == "opencl":
+        if not use_opencl():
+            print("opencl is not available")
+            return 1
+        loaders = ['opencl']
+    elif opts.engine == "dll":
+        loaders = ["dll"]
+    elif opts.engine == "cuda":
+        if not use_cuda():
+            print("cuda is not available")
+            return 1
+        loaders = ['cuda']
+    elif opts.engine == "all":
+        loaders = ['dll']
+        if use_opencl():
+            loaders.append('opencl')
+        if use_cuda():
+            loaders.append('cuda')
+    else:
+        print("unknown engine " + opts.engine)
+        return 1
+    runner = TestRunner(verbosity=opts.verbose, **test_args)
+    result = runner.run(make_suite(loaders, opts.models))
+    return 1 if result.failures or result.errors else 0
 if __name__ == "__main__":
     sys.exit(main(*sys.argv[1:]))
+    sys.exit(main())

sasmodels/modelinfo.py

-                      rc1799d3
+                      ra34b811
     other sld parameters. The volume parameters are used for calls
     to form_volume within the kernel (required for volume normalization),
     to shell_volume (for hollow shapes), and to effective_radius (for
+    to shell_volume (for hollow shapes), and to radius_effective (for
     structure factor interactions) respectively.
 …
     example, might be used to set the value of a shape parameter.
+    These values are set by :func:`make_parameter_table` and
+    Control parameters are used for variant models such as :ref:`rpa` which
+    have different cases with different parameters, as well as models
+    like :ref:`spherical_sld` with its user defined number of shells.
+    The control parameter should appear in the parameter table along with the
+    parameters it is is controlling.  For variant models, use *[CASES]* in
+    place of the parameter limits Within the parameter definition table,
+    with case names such as::
+         CASES = ["diblock copolymer", "triblock copolymer", ...]
+    This should give *limits=[[case1, case2, ...]]*, but the model loader
+    translates it to *limits=[0, len(CASES)-1]*, and adds *choices=CASES* to
+    the :class:`Parameter` definition. Note that models can use a list of
+    cases as a parameter without it being a control parameter.  Either way,
+    the parameter is sent to the model evaluator as *float(choice_num)*,
+    where choices are numbered from 0. :meth:`ModelInfo.get_hidden_parameters`
+    will determine which parameers to display.
+    The class contructor should not be called directly, but instead the
+    parameter table is built using :func:`make_parameter_table` and
     :func:`parse_parameter` therein.
     """
 …
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
     # TODO: find Fq by inspection
     info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None)
+    info.radius_effective_modes = getattr(kernel_module, 'radius_effective_modes', None)
     info.have_Fq = getattr(kernel_module, 'have_Fq', False)
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
 …
     info.source = getattr(kernel_module, 'source', [])
     info.c_code = getattr(kernel_module, 'c_code', None)
     info.effective_radius = getattr(kernel_module, 'effective_radius', None)
+    info.radius_effective = getattr(kernel_module, 'radius_effective', None)
     # TODO: check the structure of the tests
     info.tests = getattr(kernel_module, 'tests', [])
 …
     info.single = getattr(kernel_module, 'single', not callable(info.Iq))
     info.random = getattr(kernel_module, 'random', None)
-    # multiplicity info
-    control_pars = [p.id for p in parameters.kernel_parameters if p.is_control]
-    default_control = control_pars[0] if control_pars else None
-    info.control = getattr(kernel_module, 'control', default_control)
     info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore
+    # Set control flag for explicitly set parameters, e.g., in the RPA model.
+    control = getattr(kernel_module, 'control', None)
+    if control is not None:
+        parameters[control].is_control = True
     if callable(info.Iq) and parameters.has_2d:
 …
     #: *sphere*hardsphere* or *cylinder+sphere*.
     composition = None      # type: Optional[Tuple[str, List[ModelInfo]]]
-    #: Name of the control parameter for a variant model such as :ref:`rpa`.
-    #: The *control* parameter should appear in the parameter table, with
-    #: limits defined as *[CASES]*, for case names such as
-    #: *CASES = ["diblock copolymer", "triblock copolymer", ...]*.
-    #: This should give *limits=[[case1, case2, ...]]*, but the
-    #: model loader translates this to *limits=[0, len(CASES)-1]*, and adds
-    #: *choices=CASES* to the :class:`Parameter` definition. Note that
-    #: models can use a list of cases as a parameter without it being a
-    #: control parameter.  Either way, the parameter is sent to the model
-    #: evaluator as *float(choice_num)*, where choices are numbered from 0.
-    #: See also :attr:`hidden`.
-    control = None          # type: str
     #: Different variants require different parameters.  In order to show
     #: just the parameters needed for the variant selected by :attr:`control`,
 …
     #: List of options for computing the effective radius of the shape,
     #: or None if the model is not usable as a form factor model.
     effective_radius_type = None   # type: List[str]
+    radius_effective_modes = None   # type: List[str]
     #: List of C source files used to define the model.  The source files
     #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the
 …
     #: C code, either defined as a string, or in the sources.
     shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]]
+    #: Computes the effective radius of the shape given the volume parameters.
+    #: Only needed for models defined in python that can be used for
+    #: monodisperse approximation for non-dilute solutions, P@S.  The first
+    #: argument is the integer effective radius mode, with default 0.
+    radius_effective = None  # type: Union[None, Callable[[int, np.ndarray], float]]
     #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined
     #: by the parameter table.  *Iq* can be defined as a python function, or
 …
     #: will return *I(q, a, b, ...)*.  Multiplicity parameters are sent as
     #: pointers to doubles.  Constants in floating point expressions should
+    #: include the decimal point. See :mod:`generate` for more details.
+    Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    #: include the decimal point. See :mod:`generate` for more details. If
+    #: *have_Fq* is True, then Iq should return an interleaved array of
+    #: $[\sum F(q_1), \sum F^2(q_1), \ldots, \sum F(q_n), \sum F^2(q_n)]$.
+    Iq = None               # type: Union[None, str, Callable[[...], np.ndarray]]
+    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
+    Iqxy = None             # type: Union[None, str, Callable[[...], np.ndarray]]
     #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
     Iqac = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    Iqac = None             # type: Union[None, str, Callable[[...], np.ndarray]]
     #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
     Iqabc = None            # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    Iqabc = None            # type: Union[None, str, Callable[[...], np.ndarray]]
     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
     Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]]
+    Imagnetic = None        # type: Union[None, str, Callable[[...], np.ndarray]]
     #: Returns a model profile curve *x, y*.  If *profile* is defined, this
     #: curve will appear in response to the *Show* button in SasView.  Use

sasmodels/models/_spherepy.py

-                      r304c775
+                      ra34b811
 to the output of the software provided by the NIST (Kline, 2006).
 References
 ----------
 …
 John Wiley and Sons, New York, (1955)
+Source
+------
+`_spherepy.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/_spherepy.py>`_
+`sphere.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/sphere.c>`_
+Authorship and Verification
+----------------------------
+* **Author: P Kienzle**
+* **Last Modified by:**
 * **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06
+* **Source added by :** Steve King **Date:** March 25, 2019
 """
 …
 def form_volume(radius):
+    """Calculate volume for sphere"""
     return 1.333333333333333 * pi * radius ** 3
+def effective_radius(mode, radius):
+    return radius
+def radius_effective(mode, radius):
+    """Calculate R_eff for sphere"""
+    return radius if mode else 0.
 def Iq(q, sld, sld_solvent, radius):
+    """Calculate I(q) for sphere"""
     #print "q",q
     #print "sld,r",sld,sld_solvent,radius

sasmodels/models/adsorbed_layer.py

-                      r2d81cfe
+                      r0507e09
    Layers*, *Macromol. Symp.*, 190 (2002) 33-42.
+Source
+------
+`adsorbed_layer.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/adsorbed_layer.py>`_
 Authorship and Verification
 ----------------------------
 …
 * **Last Modified by:** Paul Kienzle **Date:** April 14, 2016
 * **Last Reviewed by:** Steve King **Date:** March 18, 2016
+* **Source added by :** Steve King **Date:** March 25, 2019
 """
 …
 def Iq(q, second_moment, adsorbed_amount, density_shell, radius,
        volfraction, sld_shell, sld_solvent):
+    """Return I(q) for adsorbed layer model."""
     with errstate(divide='ignore'):
         aa = ((sld_shell - sld_solvent)/density_shell * adsorbed_amount) / q
 …

SasView

Changeset bd91e8f in sasmodels

Legend:

.travis.yml

doc/guide/index.rst

doc/guide/plugin.rst

doc/guide/resolution.rst

explore/asymint.py

explore/beta/sasfit_compare.py

explore/beta/sasfit_compare_new.py

extra/build_linux.sh

extra/pylint.rc

sasmodels/bumps_model.py

sasmodels/compare.py

sasmodels/convert.py

sasmodels/core.py

sasmodels/data.py

sasmodels/details.py

sasmodels/direct_model.py

sasmodels/generate.py

sasmodels/gengauss.py

sasmodels/guyou.py

sasmodels/jitter.py

sasmodels/kernel.py

sasmodels/kernel_iq.c

sasmodels/kernelcl.py

sasmodels/kernelcuda.py

sasmodels/kerneldll.py

sasmodels/kernelpy.py

sasmodels/list_pars.py

sasmodels/mixture.py

sasmodels/model_test.py

sasmodels/modelinfo.py

sasmodels/models/_spherepy.py

sasmodels/models/adsorbed_layer.py