source: sasmodels/doc/guide/plugin.rst @ db1d9d5

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since db1d9d5 was db1d9d5, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

merge with master

  • Property mode set to 100644
File size: 51.9 KB
RevLine 
[990d8df]1.. _Writing_a_Plugin:
2
3Writing a Plugin Model
4======================
5
6Overview
7^^^^^^^^
8
9In addition to the models provided with the sasmodels package, you are free to
10create your own models.
11
12Models can be of three types:
13
14- A pure python model : Example -
15  `broadpeak.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/broad_peak.py>`_
16
17- A python model with embedded C : Example -
18  `sphere.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/sphere.py>`_
19
20- A python wrapper with separate C code : Example -
21  `cylinder.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.py>`_,
22  `cylinder.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/cylinder.c>`_
23
24When using SasView, plugin models should be saved to the SasView
25*plugin_models* folder *C:\\Users\\{username}\\.sasview\\plugin_models*
26(on Windows) or */Users/{username}/.sasview\\plugin_models* (on Mac).
27The next time SasView is started it will compile the plugin and add
28it to the list of *Plugin Models* in a FitPage.  Scripts can load
29the models from anywhere.
30
31The built-in modules are available in the *models* subdirectory
32of the sasmodels package.  For SasView on Windows, these will
33be found in *C:\\Program Files (x86)\\SasView\\sasmodels-data\\models*.
34On Mac OSX, these will be within the application bundle as
35*/Applications/SasView 4.0.app/Contents/Resources/sasmodels-data/models*.
36
37Other models are available for download from the
38`Model Marketplace <http://marketplace.sasview.org/>`_. You can contribute your
39own models to the Marketplace as well.
40
41Create New Model Files
42^^^^^^^^^^^^^^^^^^^^^^
43
44Copy the appropriate files to your plugin models directory (we recommend
45using the examples above as templates) as mymodel.py (and mymodel.c, etc)
46as required, where "mymodel" is the name for the model you are creating.
47
48*Please follow these naming rules:*
49
50- No capitalization and thus no CamelCase
51- If necessary use underscore to separate words (i.e. barbell not BarBell or
52  broad_peak not BroadPeak)
53- Do not include "model" in the name (i.e. barbell not BarBellModel)
54
55
56Edit New Model Files
57^^^^^^^^^^^^^^^^^^^^
58
59Model Contents
60..............
61
62The model interface definition is in the .py file.  This file contains:
63
64- a **model name**:
65   - this is the **name** string in the *.py* file
66   - titles should be:
67
68    - all in *lower* case
69    - without spaces (use underscores to separate words instead)
70    - without any capitalization or CamelCase
71    - without incorporating the word "model"
72    - examples: *barbell* **not** *BarBell*; *broad_peak* **not** *BroadPeak*;
73      *barbell* **not** *BarBellModel*
74
75- a **model title**:
76   - this is the **title** string in the *.py* file
77   - this is a one or two line description of the model, which will appear
78     at the start of the model documentation and as a tooltip in the SasView GUI
79
[3048ec6]80- a **short description**:
[990d8df]81   - this is the **description** string in the *.py* file
82   - this is a medium length description which appears when you click
83     *Description* on the model FitPage
84
85- a **parameter table**:
86   - this will be auto-generated from the *parameters* in the *.py* file
87
88- a **long description**:
89   - this is ReStructuredText enclosed between the r""" and """ delimiters
90     at the top of the *.py* file
91   - what you write here is abstracted into the SasView help documentation
92   - this is what other users will refer to when they want to know what
93     your model does; so please be helpful!
94
95- a **definition** of the model:
96   - as part of the **long description**
97
98- a **formula** defining the function the model calculates:
99   - as part of the **long description**
100
101- an **explanation of the parameters**:
102   - as part of the **long description**
103   - explaining how the symbols in the formula map to the model parameters
104
105- a **plot** of the function, with a **figure caption**:
106   - this is automatically generated from your default parameters
107
108- at least one **reference**:
109   - as part of the **long description**
110   - specifying where the reader can obtain more information about the model
111
112- the **name of the author**
113   - as part of the **long description**
114   - the *.py* file should also contain a comment identifying *who*
115     converted/created the model file
116
117Models that do not conform to these requirements will *never* be incorporated
118into the built-in library.
119
120
121Model Documentation
122...................
123
124The *.py* file starts with an r (for raw) and three sets of quotes
125to start the doc string and ends with a second set of three quotes.
126For example::
127
128    r"""
129    Definition
130    ----------
131
132    The 1D scattering intensity of the sphere is calculated in the following
133    way (Guinier, 1955)
134
135    .. math::
136
137        I(q) = \frac{\text{scale}}{V} \cdot \left[
138            3V(\Delta\rho) \cdot \frac{\sin(qr) - qr\cos(qr))}{(qr)^3}
139            \right]^2 + \text{background}
140
141    where *scale* is a volume fraction, $V$ is the volume of the scatterer,
142    $r$ is the radius of the sphere and *background* is the background level.
143    *sld* and *sld_solvent* are the scattering length densities (SLDs) of the
144    scatterer and the solvent respectively, whose difference is $\Delta\rho$.
145
146    You can included figures in your documentation, as in the following
147    figure for the cylinder model.
148
149    .. figure:: img/cylinder_angle_definition.jpg
150
151        Definition of the angles for oriented cylinders.
152
153    References
154    ----------
155
156    A Guinier, G Fournet, *Small-Angle Scattering of X-Rays*,
157    John Wiley and Sons, New York, (1955)
158    """
159
160This is where the FULL documentation for the model goes (to be picked up by
161the automatic documentation system).  Although it feels odd, you
162should start the documentation immediately with the **definition**---the model
163name, a brief description and the parameter table are automatically inserted
164above the definition, and the a plot of the model is automatically inserted
165before the **reference**.
166
167Figures can be included using the *figure* command, with the name
168of the *.png* file containing the figure and a caption to appear below the
169figure.  Figure numbers will be added automatically.
170
171See this `Sphinx cheat sheet <http://matplotlib.org/sampledoc/cheatsheet.html>`_
172for a quick guide to the documentation layout commands, or the
173`Sphinx Documentation <http://www.sphinx-doc.org/en/stable/>`_ for
174complete details.
175
176The model should include a **formula** written using LaTeX markup.
177The example above uses the *math* command to make a displayed equation.  You
178can also use *\$formula\$* for an inline formula. This is handy for defining
179the relationship between the model parameters and formula variables, such
180as the phrase "\$r\$ is the radius" used above.  The live demo MathJax
181page `<http://www.mathjax.org/>`_ is handy for checking that the equations
182will look like you intend.
183
184Math layout uses the `amsmath <http://www.ams.org/publications/authors/tex/amslatex>`_
185package for aligning equations (see amsldoc.pdf on that page for complete
186documentation). You will automatically be in an aligned environment, with
187blank lines separating the lines of the equation.  Place an ampersand before
188the operator on which to align.  For example::
189
190    .. math::
191
192      x + y &= 1 \\
193      y &= x - 1
194
195produces
196
197.. math::
198
199      x + y &= 1 \\
200      y &= x - 1
201
202If you need more control, use::
203
204    .. math::
205        :nowrap:
206
207
208Model Definition
209................
210
211Following the documentation string, there are a series of definitions::
212
213    name = "sphere"  # optional: defaults to the filename without .py
214
215    title = "Spheres with uniform scattering length density"
216
217    description = """\
218    P(q)=(scale/V)*[3V(sld-sld_solvent)*(sin(qr)-qr cos(qr))
219                    /(qr)^3]^2 + background
220        r: radius of sphere
221        V: The volume of the scatter
222        sld: the SLD of the sphere
223        sld_solvent: the SLD of the solvent
224    """
225
226    category = "shape:sphere"
227
228    single = True   # optional: defaults to True
229
230    opencl = False  # optional: defaults to False
231
232    structure_factor = False  # optional: defaults to False
233
234**name = "mymodel"** defines the name of the model that is shown to the user.
[3048ec6]235If it is not provided it will use the name of the model file. The name must
236be a valid variable name, starting with a letter and contains only letters,
237numbers or underscore.  Spaces, dashes, and other symbols are not permitted.
[990d8df]238
239**title = "short description"** is short description of the model which
240is included after the model name in the automatically generated documentation.
241The title can also be used for a tooltip.
242
243**description = """doc string"""** is a longer description of the model. It
244shows up when you press the "Description" button of the SasView FitPage.
245It should give a brief description of the equation and the parameters
246without the need to read the entire model documentation. The triple quotes
247allow you to write the description over multiple lines. Keep the lines
248short since the GUI will wrap each one separately if they are too long.
249**Make sure the parameter names in the description match the model definition!**
250
251**category = "shape:sphere"** defines where the model will appear in the
252model documentation.  In this example, the model will appear alphabetically
253in the list of spheroid models in the *Shape* category.
254
255**single = True** indicates that the model can be run using single
256precision floating point values.  Set it to False if the numerical
257calculation for the model is unstable, which is the case for about 20 of
258the built in models.  It is worthwhile modifying the calculation to support
259single precision, allowing models to run up to 10 times faster.  The
260section `Test_Your_New_Model`_  describes how to compare model values for
261single vs. double precision so you can decide if you need to set
262single to False.
263
264**opencl = False** indicates that the model should not be run using OpenCL.
265This may be because the model definition includes code that cannot be
266compiled for the GPU (for example, goto statements).  It can also be used
267for large models which can't run on most GPUs.  This flag has not been
268used on any of the built in models; models which were failing were
269streamlined so this flag was not necessary.
270
271**structure_factor = True** indicates that the model can be used as a
272structure factor to account for interactions between particles.  See
273`Form_Factors`_ for more details.
274
[9d8a027]275**model_info = ...** lets you define a model directly, for example, by
276loading and modifying existing models.  This is done implicitly by
277:func:`sasmodels.core.load_model_info`, which can create a mixture model
278from a pair of existing models.  For example::
279
280    from sasmodels.core import load_model_info
281    model_info = load_model_info('sphere+cylinder')
282
283See :class:`sasmodels.modelinfo.ModelInfo` for details about the model
284attributes that are defined.
285
[990d8df]286Model Parameters
287................
288
289Next comes the parameter table.  For example::
290
291    # pylint: disable=bad-whitespace, line-too-long
292    #   ["name",        "units", default, [min, max], "type",    "description"],
293    parameters = [
294        ["sld",         "1e-6/Ang^2",  1, [-inf, inf], "sld",    "Layer scattering length density"],
295        ["sld_solvent", "1e-6/Ang^2",  6, [-inf, inf], "sld",    "Solvent scattering length density"],
296        ["radius",      "Ang",        50, [0, inf],    "volume", "Sphere radius"],
297    ]
298    # pylint: enable=bad-whitespace, line-too-long
299
300**parameters = [["name", "units", default, [min,max], "type", "tooltip"],...]**
301defines the parameters that form the model.
302
303**Note: The order of the parameters in the definition will be the order of the
[31fc4ad]304parameters in the user interface and the order of the parameters in Fq(), Iq(),
[e5cb3df]305Iqac(), Iqabc(), radius_effective(), form_volume() and shell_volume().
[31fc4ad]306And** *scale* **and** *background* **parameters are implicit to all models,
307so they do not need to be included in the parameter table.**
[990d8df]308
309- **"name"** is the name of the parameter shown on the FitPage.
310
[3048ec6]311  - the name must be a valid variable name, starting with a letter and
312    containing only letters, numbers and underscore.
313
[990d8df]314  - parameter names should follow the mathematical convention; e.g.,
315    *radius_core* not *core_radius*, or *sld_solvent* not *solvent_sld*.
316
317  - model parameter names should be consistent between different models,
318    so *sld_solvent*, for example, should have exactly the same name
319    in every model.
320
321  - to see all the parameter names currently in use, type the following in the
322    python shell/editor under the Tools menu::
323
324       import sasmodels.list_pars
325       sasmodels.list_pars.list_pars()
326
327    *re-use* as many as possible!!!
328
329  - use "name[n]" for multiplicity parameters, where *n* is the name of
330    the parameter defining the number of shells/layers/segments, etc.
331
332- **"units"** are displayed along with the parameter name
333
334  - every parameter should have units; use "None" if there are no units.
335
336  - **sld's should be given in units of 1e-6/Ang^2, and not simply
337    1/Ang^2 to be consistent with the builtin models.  Adjust your formulas
338    appropriately.**
339
340  - fancy units markup is available for some units, including::
341
342        Ang, 1/Ang, 1/Ang^2, 1e-6/Ang^2, degrees, 1/cm, Ang/cm, g/cm^3, mg/m^2
343
344  - the list of units is defined in the variable *RST_UNITS* within
345    `sasmodels/generate.py <https://github.com/SasView/sasmodels/tree/master/sasmodels/generate.py>`_
346
347    - new units can be added using the macros defined in *doc/rst_prolog*
348      in the sasmodels source.
349    - units should be properly formatted using sub-/super-scripts
350      and using negative exponents instead of the / operator, though
351      the unit name should use the / operator for consistency.
352    - please post a message to the SasView developers mailing list with your changes.
353
354- **default** is the initial value for the parameter.
355
356  - **the parameter default values are used to auto-generate a plot of
357    the model function in the documentation.**
358
359- **[min, max]** are the lower and upper limits on the parameter.
360
361  - lower and upper limits can be any number, or *-inf* or *inf*.
362
363  - the limits will show up as the default limits for the fit making it easy,
364    for example, to force the radius to always be greater than zero.
365
366  - these are hard limits defining the valid range of parameter values;
367    polydisperity distributions will be truncated at the limits.
368
369- **"type"** can be one of: "", "sld", "volume", or "orientation".
370
371  - "sld" parameters can have magnetic moments when fitting magnetic models;
372    depending on the spin polarization of the beam and the $q$ value being
373    examined, the effective sld for that material will be used to compute the
374    scattered intensity.
375
[31fc4ad]376  - "volume" parameters are passed to Fq(), Iq(), Iqac(), Iqabc(), form_volume()
377    and shell_volume(), and have polydispersity loops generated automatically.
[990d8df]378
[108e70e]379  - "orientation" parameters are not passed, but instead are combined with
380    orientation dispersity to translate *qx* and *qy* to *qa*, *qb* and *qc*.
381    These parameters should appear at the end of the table with the specific
382    names *theta*, *phi* and for asymmetric shapes *psi*, in that order.
[990d8df]383
[9844c3a]384Some models will have integer parameters, such as number of pearls in the
385pearl necklace model, or number of shells in the multi-layer vesicle model.
386The optimizers in BUMPS treat all parameters as floating point numbers which
387can take arbitrary values, even for integer parameters, so your model should
388round the incoming parameter value to the nearest integer inside your model
[e5cb3df]389you should round to the nearest integer.  In C code, you can do this using:
390
391.. code-block:: c
[9844c3a]392
393    static double
394    Iq(double q, ..., double fp_n, ...)
395    {
396        int n = (int)(fp_n + 0.5);
397        ...
398    }
399
400in python::
401
402    def Iq(q, ..., n, ...):
403        n = int(n+0.5)
404        ...
405
[3048ec6]406Derivative based optimizers such as Levenberg-Marquardt will not work
[9844c3a]407for integer parameters since the partial derivative is always zero, but
408the remaining optimizers (DREAM, differential evolution, Nelder-Mead simplex)
409will still function.
410
[990d8df]411Model Computation
412.................
413
414Models can be defined as pure python models, or they can be a mixture of
415python and C models.  C models are run on the GPU if it is available,
416otherwise they are compiled and run on the CPU.
417
418Models are defined by the scattering kernel, which takes a set of parameter
419values defining the shape, orientation and material, and returns the
420expected scattering. Polydispersity and angular dispersion are defined
421by the computational infrastructure.  Any parameters defined as "volume"
422parameters are polydisperse, with polydispersity defined in proportion
423to their value.  "orientation" parameters use angular dispersion defined
424in degrees, and are not relative to the current angle.
425
426Based on a weighting function $G(x)$ and a number of points $n$, the
427computed value is
428
429.. math::
430
431     \hat I(q)
432     = \frac{\int G(x) I(q, x)\,dx}{\int G(x) V(x)\,dx}
433     \approx \frac{\sum_{i=1}^n G(x_i) I(q,x_i)}{\sum_{i=1}^n G(x_i) V(x_i)}
434
[3048ec6]435That is, the individual models do not need to include polydispersity
[990d8df]436calculations, but instead rely on numerical integration to compute the
[108e70e]437appropriately smeared pattern.
[990d8df]438
[2015f02]439Each .py file also contains a function::
440
441        def random():
442        ...
[31fc4ad]443
444This function provides a model-specific random parameter set which shows model
445features in the USANS to SANS range.  For example, core-shell sphere sets the
446outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q
447value for the transition from flat to falling.  It then uses a beta distribution
448to set the percentage of the shape which is shell, giving a preference for very
449thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp
450should show a reasonable variety of curves over the default sascomp q range.
451The parameter set is returned as a dictionary of `{parameter: value, ...}`.
452Any model parameters not included in the dictionary will default according to
[2015f02]453the code in the `_randomize_one()` function from sasmodels/compare.py.
454
[990d8df]455Python Models
456.............
457
458For pure python models, define the *Iq* function::
459
460      import numpy as np
461      from numpy import cos, sin, ...
462
463      def Iq(q, par1, par2, ...):
464          return I(q, par1, par2, ...)
465      Iq.vectorized = True
466
467The parameters *par1, par2, ...* are the list of non-orientation parameters
468to the model in the order that they appear in the parameter table.
[3048ec6]469**Note that the auto-generated model file uses** *x* **rather than** *q*.
[990d8df]470
471The *.py* file should import trigonometric and exponential functions from
472numpy rather than from math.  This lets us evaluate the model for the whole
473range of $q$ values at once rather than looping over each $q$ separately in
474python.  With $q$ as a vector, you cannot use if statements, but must instead
475do tricks like
476
477::
478
479     a = x*q*(q>0) + y*q*(q<=0)
480
481or
482
483::
484
485     a = np.empty_like(q)
486     index = q>0
487     a[index] = x*q[index]
488     a[~index] = y*q[~index]
489
490which sets $a$ to $q \cdot x$ if $q$ is positive or $q \cdot y$ if $q$
491is zero or negative. If you have not converted your function to use $q$
492vectors, you can set the following and it will only receive one $q$
493value at a time::
494
495    Iq.vectorized = False
496
497Return np.NaN if the parameters are not valid (e.g., cap_radius < radius in
498barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be
499ignored, and not included in the calculation of the weighted polydispersity.
500
501Models should define *form_volume(par1, par2, ...)* where the parameter
502list includes the *volume* parameters in order.  This is used for a weighted
[e5cb3df]503volume normalization so that scattering is on an absolute scale.  For
504solid shapes, the *I(q)* function should use *form_volume* squared
505as its scale factor.  If *form_volume* is not defined, then the default
506*form_volume = 1.0* will be used.
[990d8df]507
[31fc4ad]508Hollow shapes, where the volume fraction of particle corresponds to the
509material in the shell rather than the volume enclosed by the shape, must
510also define a *shell_volume(par1, par2, ...)* function.  The parameters
[e5cb3df]511are the same as for *form_volume*.  Here the *I(q)* function should use
512*shell_volume* squared instead of *form_volume* squared so that the scale
513parameter corresponds to the volume fraction of material in the sample.
514The structure factor calculation needs the volume fraction of the filled
515shapes for its calculation, so the volume fraction parameter in the model
516is automatically scaled by *form_volume/shell_volume* prior to calling the
517structure factor.
[31fc4ad]518
[2fe39d1]519**Note: Pure python models do not yet support direct computation of the**
520**average of $F(q)$ and $F^2(q)$. Neither do they support orientational**
521**distributions or magnetism (use C models if these are required).**
[31fc4ad]522
[990d8df]523Embedded C Models
524.................
525
526Like pure python models, inline C models need to define an *Iq* function::
527
528    Iq = """
529        return I(q, par1, par2, ...);
530    """
531
[e5cb3df]532This expands into the equivalent C code:
533
534.. code-block:: c
[990d8df]535
536    double Iq(double q, double par1, double par2, ...);
537    double Iq(double q, double par1, double par2, ...)
538    {
539        return I(q, par1, par2, ...);
540    }
541
542*form_volume* defines the volume of the shape. As in python models, it
543includes only the volume parameters.
544
[e5cb3df]545*shell_volume* defines the volume of the shell for hollow shapes. As in
[31fc4ad]546python models, it includes only the volume parameters.
547
[990d8df]548**source=['fn.c', ...]** includes the listed C source files in the
[108e70e]549program before *Iq* and *form_volume* are defined. This allows you to
[ef85a09]550extend the library of C functions available to your model.
551
552*c_code* includes arbitrary C code into your kernel, which can be
553handy for defining helper functions for *Iq* and *form_volume*. Note that
[108e70e]554you can put the full function definition for *Iq* and *form_volume*
[ef85a09]555(include function declaration) into *c_code* as well, or put them into an
556external C file and add that file to the list of sources.
[990d8df]557
558Models are defined using double precision declarations for the
559parameters and return values.  When a model is run using single
560precision or long double precision, each variable is converted
561to the target type, depending on the precision requested.
562
563**Floating point constants must include the decimal point.**  This allows us
564to convert values such as 1.0 (double precision) to 1.0f (single precision)
565so that expressions that use these values are not promoted to double precision
566expressions.  Some graphics card drivers are confused when functions
567that expect floating point values are passed integers, such as 4*atan(1); it
568is safest to not use integers in floating point expressions.  Even better,
569use the builtin constant M_PI rather than 4*atan(1); it is faster and smaller!
570
571The C model operates on a single $q$ value at a time.  The code will be
572run in parallel across different $q$ values, either on the graphics card
573or the processor.
574
575Rather than returning NAN from Iq, you must define the *INVALID(v)*.  The
576*v* parameter lets you access all the parameters in the model using
[e5cb3df]577*v.par1*, *v.par2*, etc. For example:
578
579.. code-block:: c
[990d8df]580
581    #define INVALID(v) (v.bell_radius < v.radius)
582
[ef85a09]583The INVALID define can go into *Iq*, or *c_code*, or an external C file
584listed in *source*.
585
[31fc4ad]586Structure Factors
587.................
588
589Structure factor calculations may need the underlying $<F(q)>$ and $<F^2(q)>$
590rather than $I(q)$.  This is used to compute $\beta = <F(q)>^2/<F^2(q)>$ in
591the decoupling approximation to the structure factor.
592
593Instead of defining the *Iq* function, models can define *Fq* as
[e5cb3df]594something like:
595
596.. code-block:: c
[31fc4ad]597
598    double Fq(double q, double *F1, double *F2, double par1, double par2, ...);
599    double Fq(double q, double *F1, double *F2, double par1, double par2, ...)
600    {
601        // Polar integration loop over all orientations.
602        ...
603        *F1 = 1e-2 * total_F1 * contrast * volume;
604        *F2 = 1e-4 * total_F2 * square(contrast * volume);
605        return I(q, par1, par2, ...);
606    }
607
608If the volume fraction scale factor is built into the model (as occurs for
609the vesicle model, for example), then scale *F1* by $\surd V_f$ so that
610$\beta$ is computed correctly.
611
612Structure factor calculations are not yet supported for oriented shapes.
613
614Note: only available as a separate C file listed in *source*, or within
615a *c_code* block within the python model definition file.
616
[108e70e]617Oriented Shapes
618...............
619
620If the scattering is dependent on the orientation of the shape, then you
621will need to include *orientation* parameters *theta*, *phi* and *psi*
[7e6bc45e]622at the end of the parameter table.  As described in the section
623:ref:`orientation`, the individual $(q_x, q_y)$ points on the detector will
624be rotated into $(q_a, q_b, q_c)$ points relative to the sample in its
625canonical orientation with $a$-$b$-$c$ aligned with $x$-$y$-$z$ in the
626laboratory frame and beam travelling along $-z$.
627
[e5cb3df]628The oriented C model (oriented pure Python models are not supported)
[2fe39d1]629is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where
[108e70e]630*par1*, etc. are the parameters to the model.  If the shape is rotationally
631symmetric about *c* then *psi* is not needed, and the model is called
632as *Iqac(qab, qc, par1, par2, ...)*.  In either case, the orientation
633parameters are not included in the function call.
634
635For 1D oriented shapes, an integral over all angles is usually needed for
[b85227d]636the *Iq* function. Given symmetry and the substitution $u = \cos(\alpha)$,
[108e70e]637$du = -\sin(\alpha)\,d\alpha$ this becomes
638
639.. math::
640
[b85227d]641    I(q) &= \frac{1}{4\pi} \int_{-\pi/2}^{pi/2} \int_{-pi}^{pi}
642            F(q_a, q_b, q_c)^2 \sin(\alpha)\,d\beta\,d\alpha \\
643        &= \frac{8}{4\pi} \int_{0}^{pi/2} \int_{0}^{\pi/2}
644            F^2 \sin(\alpha)\,d\beta\,d\alpha \\
645        &= \frac{8}{4\pi} \int_1^0 \int_{0}^{\pi/2} - F^2 \,d\beta\,du \\
646        &= \frac{8}{4\pi} \int_0^1 \int_{0}^{\pi/2} F^2 \,d\beta\,du
647
648for
649
650.. math::
651
652    q_a &= q \sin(\alpha)\sin(\beta) = q \sqrt{1-u^2} \sin(\beta) \\
653    q_b &= q \sin(\alpha)\cos(\beta) = q \sqrt{1-u^2} \cos(\beta) \\
654    q_c &= q \cos(\alpha) = q u
[108e70e]655
656Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the
[e5cb3df]657numerical integration is then:
658
659.. code-block:: c
[108e70e]660
661    double outer_sum = 0.0;
662    for (int i = 0; i < GAUSS_N; i++) {
663        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5;
664        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
665        const double qc = cos_alpha * q;
666        double inner_sum = 0.0;
667        for (int j = 0; j < GAUSS_N; j++) {
668            const double beta = M_PI_4 * GAUSS_Z[j] + M_PI_4;
669            double sin_beta, cos_beta;
670            SINCOS(beta, sin_beta, cos_beta);
671            const double qa = sin_alpha * sin_beta * q;
[b85227d]672            const double qb = sin_alpha * cos_beta * q;
673            const double form = Fq(qa, qb, qc, ...);
674            inner_sum += GAUSS_W[j] * form * form;
[108e70e]675        }
676        outer_sum += GAUSS_W[i] * inner_sum;
677    }
678    outer_sum *= 0.25; // = 8/(4 pi) * outer_sum * (pi/2) / 4
679
680The *z* values for the Gauss-Legendre integration extends from -1 to 1, so
681the double sum of *w[i]w[j]* explains the factor of 4.  Correcting for the
682average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$.  The $8/(4 \pi)$
683factor comes from the integral over the quadrant.  With less symmetry (eg.,
684in the bcc and fcc paracrystal models), then an integral over the entire
685sphere may be necessary.
686
687For simpler models which are rotationally symmetric a single integral
688suffices:
689
690.. math::
691
[b85227d]692    I(q) &= \frac{1}{\pi}\int_{-\pi/2}^{\pi/2}
693            F(q_{ab}, q_c)^2 \sin(\alpha)\,d\alpha/\pi \\
694        &= \frac{2}{\pi} \int_0^1 F^2\,du
695
696for
697
698.. math::
699
700    q_{ab} &= q \sin(\alpha) = q \sqrt{1 - u^2} \\
701    q_c &= q \cos(\alpha) = q u
702
[108e70e]703
704with integration loop::
705
706    double sum = 0.0;
707    for (int i = 0; i < GAUSS_N; i++) {
708        const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5;
709        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
710        const double qab = sin_alpha * q;
[b85227d]711        const double qc = cos_alpha * q;
712        const double form = Fq(qab, qc, ...);
713        sum += GAUSS_W[j] * form * form;
[108e70e]714    }
715    sum *= 0.5; // = 2/pi * sum * (pi/2) / 2
716
717Magnetism
718.........
719
720Magnetism is supported automatically for all shapes by modifying the
721effective SLD of particle according to the Halpern-Johnson vector
[c654160]722describing the interaction between neutron spin and magnetic field.  All
[108e70e]723parameters marked as type *sld* in the parameter table are treated as
724possibly magnetic particles with magnitude *M0* and direction
725*mtheta* and *mphi*.  Polarization parameters are also provided
726automatically for magnetic models to set the spin state of the measurement.
727
728For more complicated systems where magnetism is not uniform throughout
729the individual particles, you will need to write your own models.
730You should not mark the nuclear sld as type *sld*, but instead leave
731them unmarked and provide your own magnetism and polarization parameters.
732For 2D measurements you will need $(q_x, q_y)$ values for the measurement
733to compute the proper magnetism and orientation, which you can implement
734using *Iqxy(qx, qy, par1, par2, ...)*.
735
[2fe39d1]736**Note: Magnetism is not supported in pure Python models.**
737
[990d8df]738Special Functions
739.................
740
741The C code follows the C99 standard, with the usual math functions,
742as defined in
743`OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_.
744This includes the following:
745
746    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E:
747        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$
[d0dc9a3]748    exp, log, pow(x,y), expm1, log1p, sqrt, cbrt:
749        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$,
750        $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x)
751        are accurate across all $x$, including $x$ very close to zero.
[990d8df]752    sin, cos, tan, asin, acos, atan:
753        Trigonometry functions and inverses, operating on radians.
754    sinh, cosh, tanh, asinh, acosh, atanh:
755        Hyperbolic trigonometry functions.
756    atan2(y,x):
757        Angle from the $x$\ -axis to the point $(x,y)$, which is equal to
758        $\tan^{-1}(y/x)$ corrected for quadrant.  That is, if $x$ and $y$ are
759        both negative, then atan2(y,x) returns a value in quadrant III where
760        atan(y/x) would return a value in quadrant I. Similarly for
761        quadrants II and IV when $x$ and $y$ have opposite sign.
[d0dc9a3]762    fabs(x), fmin(x,y), fmax(x,y), trunc, rint:
[990d8df]763        Floating point functions.  rint(x) returns the nearest integer.
764    NAN:
765        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that
766        you cannot use :code:`x == NAN` to test for NaN values since that
[d0dc9a3]767        will always return false.  NAN does not equal NAN!  The alternative,
768        :code:`x != x` may fail if the compiler optimizes the test away.
[990d8df]769    INFINITY:
770        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x)
771        to test for finite and not NaN.
772    erf, erfc, tgamma, lgamma:  **do not use**
773        Special functions that should be part of the standard, but are missing
[fba9ca0]774        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma
775        and sas_lgamma instead (see below).
[990d8df]776
777Some non-standard constants and functions are also provided:
778
779    M_PI_180, M_4PI_3:
780        $\frac{\pi}{180}$, $\frac{4\pi}{3}$
781    SINCOS(x, s, c):
782        Macro which sets s=sin(x) and c=cos(x). The variables *c* and *s*
783        must be declared first.
784    square(x):
785        $x^2$
786    cube(x):
787        $x^3$
788    sas_sinx_x(x):
789        $\sin(x)/x$, with limit $\sin(0)/0 = 1$.
790    powr(x, y):
791        $x^y$ for $x \ge 0$; this is faster than general $x^y$ on some GPUs.
792    pown(x, n):
793        $x^n$ for $n$ integer; this is faster than general $x^n$ on some GPUs.
794    FLOAT_SIZE:
795        The number of bytes in a floating point value.  Even though all
796        variables are declared double, they may be converted to single
797        precision float before running. If your algorithm depends on
798        precision (which is not uncommon for numerical algorithms), use
799        the following::
800
801            #if FLOAT_SIZE>4
802            ... code for double precision ...
803            #else
804            ... code for single precision ...
805            #endif
806    SAS_DOUBLE:
807        A replacement for :code:`double` so that the declared variable will
808        stay double precision; this should generally not be used since some
809        graphics cards do not support double precision.  There is no provision
810        for forcing a constant to stay double precision.
811
812The following special functions and scattering calculations are defined in
813`sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_.
814These functions have been tuned to be fast and numerically stable down
815to $q=0$ even in single precision.  In some cases they work around bugs
816which appear on some platforms but not others, so use them where needed.
817Add the files listed in :code:`source = ["lib/file.c", ...]` to your *model.py*
818file in the order given, otherwise these functions will not be available.
819
820    polevl(x, c, n):
821        Polynomial evaluation $p(x) = \sum_{i=0}^n c_i x^i$ using Horner's
822        method so it is faster and more accurate.
823
824        $c = \{c_n, c_{n-1}, \ldots, c_0 \}$ is the table of coefficients,
825        sorted from highest to lowest.
826
827        :code:`source = ["lib/polevl.c", ...]` (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
828
829    p1evl(x, c, n):
830        Evaluation of normalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$
831        using Horner's method so it is faster and more accurate.
832
833        $c = \{c_{n-1}, c_{n-2} \ldots, c_0 \}$ is the table of coefficients,
834        sorted from highest to lowest.
835
836        :code:`source = ["lib/polevl.c", ...]`
[870a2f4]837        (`polevl.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
[990d8df]838
839    sas_gamma(x):
[30b60d2]840        Gamma function sas_gamma\ $(x) = \Gamma(x)$.
[990d8df]841
[fba9ca0]842        The standard math function, tgamma(x), is unstable for $x < 1$
[990d8df]843        on some platforms.
844
[870a2f4]845        :code:`source = ["lib/sas_gamma.c", ...]`
846        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_)
[990d8df]847
[fba9ca0]848    sas_gammaln(x):
849        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$.
850
851        The standard math function, lgamma(x), is incorrect for single
852        precision on some platforms.
853
854        :code:`source = ["lib/sas_gammainc.c", ...]`
855        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
856
857    sas_gammainc(a, x), sas_gammaincc(a, x):
858        Incomplete gamma function
859        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$
860        and complementary incomplete gamma function
861        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$
862
863        :code:`source = ["lib/sas_gammainc.c", ...]`
864        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
865
[990d8df]866    sas_erf(x), sas_erfc(x):
867        Error function
[30b60d2]868        sas_erf\ $(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
[990d8df]869        and complementary error function
[30b60d2]870        sas_erfc\ $(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
[990d8df]871
872        The standard math functions erf(x) and erfc(x) are slower and broken
873        on some platforms.
874
875        :code:`source = ["lib/polevl.c", "lib/sas_erf.c", ...]`
[870a2f4]876        (`sas_erf.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_erf.c>`_)
[990d8df]877
878    sas_J0(x):
[30b60d2]879        Bessel function of the first kind sas_J0\ $(x)=J_0(x)$ where
[990d8df]880        $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$.
881
882        The standard math function j0(x) is not available on all platforms.
883
884        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", ...]`
[870a2f4]885        (`sas_J0.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J0.c>`_)
[990d8df]886
887    sas_J1(x):
[30b60d2]888        Bessel function of the first kind  sas_J1\ $(x)=J_1(x)$ where
[990d8df]889        $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$.
890
891        The standard math function j1(x) is not available on all platforms.
892
893        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
[870a2f4]894        (`sas_J1.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
[990d8df]895
896    sas_JN(n, x):
[30b60d2]897        Bessel function of the first kind and integer order $n$,
898        sas_JN\ $(n, x) =J_n(x)$ where
[990d8df]899        $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$.
[30b60d2]900        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively.
[990d8df]901
[57c609b]902        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100].
903
[990d8df]904        The standard math function jn(n, x) is not available on all platforms.
905
906        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c", ...]`
[870a2f4]907        (`sas_JN.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_JN.c>`_)
[990d8df]908
909    sas_Si(x):
[30b60d2]910        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
[990d8df]911
[57c609b]912        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100].
913
[990d8df]914        This function uses Taylor series for small and large arguments:
915
[57c609b]916        For large arguments use the following Taylor series,
[990d8df]917
918        .. math::
919
920             \text{Si}(x) \sim \frac{\pi}{2}
921             - \frac{\cos(x)}{x}\left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right)
922             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)
923
[94bfa42]924        For small arguments,
[990d8df]925
926        .. math::
927
928           \text{Si}(x) \sim x
929           - \frac{x^3}{3\times 3!} + \frac{x^5}{5 \times 5!} - \frac{x^7}{7 \times 7!}
930           + \frac{x^9}{9\times 9!} - \frac{x^{11}}{11\times 11!}
931
932        :code:`source = ["lib/Si.c", ...]`
[f796469]933        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_Si.c>`_)
[990d8df]934
935    sas_3j1x_x(x):
936        Spherical Bessel form
[30b60d2]937        sph_j1c\ $(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
[990d8df]938        with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical
939        Bessel function of the first kind and first order.
940
941        This function uses a Taylor series for small $x$ for numerical accuracy.
942
943        :code:`source = ["lib/sas_3j1x_x.c", ...]`
[870a2f4]944        (`sas_3j1x_x.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_3j1x_x.c>`_)
[990d8df]945
946
947    sas_2J1x_x(x):
[30b60d2]948        Bessel form sas_J1c\ $(x) = 2 J_1(x)/x$, with a limiting value
[990d8df]949        of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind
950        and first order.
951
952        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
[870a2f4]953        (`sas_J1.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
[990d8df]954
955
956    Gauss76Z[i], Gauss76Wt[i]:
957        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively,
958        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$.
959
960        Similar arrays are available in :code:`gauss20.c` for 20-point
961        quadrature and in :code:`gauss150.c` for 150-point quadrature.
[d0dc9a3]962        The macros :code:`GAUSS_N`, :code:`GAUSS_Z` and :code:`GAUSS_W` are
963        defined so that you can change the order of the integration by
964        selecting an different source without touching the C code.
[990d8df]965
966        :code:`source = ["lib/gauss76.c", ...]`
[870a2f4]967        (`gauss76.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/gauss76.c>`_)
[990d8df]968
969
970
971Problems with C models
972......................
973
974The graphics processor (GPU) in your computer is a specialized computer tuned
975for certain kinds of problems.  This leads to strange restrictions that you
976need to be aware of.  Your code may work fine on some platforms or for some
977models, but then return bad values on other platforms.  Some examples of
978particular problems:
979
980  **(1) Code is too complex, or uses too much memory.** GPU devices only
981  have a limited amount of memory available for each processor. If you run
982  programs which take too much memory, then rather than running multiple
983  values in parallel as it usually does, the GPU may only run a single
984  version of the code at a time, making it slower than running on the CPU.
985  It may fail to run on some platforms, or worse, cause the screen to go
986  blank or the system to reboot.
987
988  **(2) Code takes too long.** Because GPU devices are used for the computer
989  display, the OpenCL drivers are very careful about the amount of time they
990  will allow any code to run. For example, on OS X, the model will stop
991  running after 5 seconds regardless of whether the computation is complete.
992  You may end up with only some of your 2D array defined, with the rest
993  containing random data. Or it may cause the screen to go blank or the
994  system to reboot.
995
996  **(3) Memory is not aligned**. The GPU hardware is specialized to operate
997  on multiple values simultaneously. To keep the GPU simple the values in
998  memory must be aligned with the different GPU compute engines. Not
999  following these rules can lead to unexpected values being loaded into
1000  memory, and wrong answers computed. The conclusion from a very long and
1001  strange debugging session was that any arrays that you declare in your
[e5cb3df]1002  model should be a multiple of four. For example:
1003
1004  .. code-block:: c
[990d8df]1005
1006      double Iq(q, p1, p2, ...)
1007      {
1008          double vector[8];  // Only going to use seven slots, but declare 8
1009          ...
1010      }
1011
1012The first step when your model is behaving strangely is to set
1013**single=False**. This automatically restricts the model to only run on the
1014CPU, or on high-end GPU cards. There can still be problems even on high-end
1015cards, so you can force the model off the GPU by setting **opencl=False**.
1016This runs the model as a normal C program without any GPU restrictions so
1017you know that strange results are probably from your code rather than the
1018environment. Once the code is debugged, you can compare your output to the
1019output on the GPU.
1020
1021Although it can be difficult to get your model to work on the GPU, the reward
1022can be a model that runs 1000x faster on a good card.  Even your laptop may
1023show a 50x improvement or more over the equivalent pure python model.
1024
1025
1026.. _Form_Factors:
1027
1028Form Factors
1029............
1030
1031Away from the dilute limit you can estimate scattering including
1032particle-particle interactions using $I(q) = P(q)*S(q)$ where $P(q)$
1033is the form factor and $S(q)$ is the structure factor.  The simplest
1034structure factor is the *hardsphere* interaction, which
1035uses the effective radius of the form factor as an input to the structure
[e5cb3df]1036factor model.  The effective radius is the weighted average over all
1037values of the shape in polydisperse systems.
1038
1039There can be many notions of effective radius, depending on the shape.  For
1040a sphere it is clearly just the radius, but for an ellipsoid of revolution
1041we provide average curvature, equivalent sphere radius, minimum radius and
1042maximum radius.  These options are listed as *radius_effective_modes* in
1043the python model defintion, and must be computed by the *radius_effective*
1044function which takes the *radius_effective_mode* parameter as an integer,
1045along with the various model parameters.  Unlike normal C/Python arrays,
1046the first mode is 1, the second is 2, etc.  Mode 0 indicates that the
1047effective radius from the shape is to be ignored in favour of the the
1048effective radius parameter in the structure factor model.
1049
1050
1051Consider the core-shell sphere, which defines the following effective radius
1052modes in the python model::
1053
1054    radius_effective_modes = [
1055        "outer radius",
1056        "core radius",
1057    ]
[990d8df]1058
[e5cb3df]1059and the following function in the C-file for the model:
[990d8df]1060
[e5cb3df]1061.. code-block:: c
[990d8df]1062
[e5cb3df]1063    static double
1064    radius_effective(int mode, double radius, double thickness)
1065    {
1066        switch (mode) {
1067            case 0: return radius + thickness;
1068            case 1: return radius;
1069            default: return 0.;
1070        }
1071    }
1072
1073    static double
1074    form_volume(double radius, double thickness)
1075    {
1076        return M_4PI_3 * cube(radius + thickness);
1077    }
[990d8df]1078
[e5cb3df]1079Given polydispersity over *(r1, r2, ..., rm)* in radius and *(t1, t2, ..., tn)*
1080in thickness, *radius_effective* is called over a mesh grid covering all
1081possible combinations of radius and thickness, with a single *(ri, tj)* pair
1082in each call. The weights each of these results according to the
1083polydispersity distributions and calls the structure factor with the average
1084effective radius.  Similarly, for *form_volume*.
[990d8df]1085
[e5cb3df]1086Hollow models have an additional volume ratio which is needed to scale the
1087structure factor.  The structure factor uses the volume fraction of the filled
1088particles as part of its density estimate, but the scale factor for the
1089scattering intensity (as non-solvent volume fraction / volume) is determined
1090by the shell volume only.  Therefore the *shell_volume* function is
1091needed to compute the form:shell volume ratio, which then scales the
1092*volfraction* parameter prior to calling the structure factor calculator.
1093In the case of a hollow sphere, this would be:
1094
1095.. code-block:: c
1096
1097    static double
1098    shell_volume(double radius, double thickness)
1099    {
1100        double whole = M_4PI_3 * cube(radius + thickness);
1101        double core = M_4PI_3 * cube(radius);
1102        return whole - core;
1103    }
[990d8df]1104
[e5cb3df]1105If *shell_volume* is not present, then *form_volume* and *shell_volume* are
1106assumed to be equal, and the shape is considered solid.
[990d8df]1107
1108Unit Tests
1109..........
1110
1111THESE ARE VERY IMPORTANT. Include at least one test for each model and
1112PLEASE make sure that the answer value is correct (i.e. not a random number).
1113
1114::
1115
1116    tests = [
1117        [{}, 0.2, 0.726362],
1118        [{"scale": 1., "background": 0., "sld": 6., "sld_solvent": 1.,
1119          "radius": 120., "radius_pd": 0.2, "radius_pd_n":45},
1120         0.2, 0.228843],
[304c775]1121        [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45},
1122         0.1, None, None, 120., None, 1.],  # q, F, F^2, R_eff, V, form:shell
[81751c2]1123        [{"@S": "hardsphere"}, 0.1, None],
[990d8df]1124    ]
1125
1126
[304c775]1127**tests=[[{parameters}, q, Iq], ...]** is a list of lists.
[990d8df]1128Each list is one test and contains, in order:
1129
1130- a dictionary of parameter values. This can be *{}* using the default
1131  parameters, or filled with some parameters that will be different from the
1132  default, such as *{"radius":10.0, "sld":4}*. Unlisted parameters will
1133  be given the default values.
1134- the input $q$ value or tuple of $(q_x, q_y)$ values.
1135- the output $I(q)$ or $I(q_x,q_y)$ expected of the model for the parameters
1136  and input value given.
1137- input and output values can themselves be lists if you have several
1138  $q$ values to test for the same model parameters.
[304c775]1139- for testing effective radius, volume and form:shell volume ratio, use the
1140  extended form of the tests results, with *None, None, R_eff, V, V_r*
1141  instead of *Iq*.  This calls the kernel *Fq* function instead of *Iq*.
1142- for testing F and F^2 (used for beta approximation) do the same as the
1143  effective radius test, but include values for the first two elements,
1144  $<F(q)>$ and $<F^2(q)>$.
[81751c2]1145- for testing interaction between form factor and structure factor, specify
1146  the structure factor name in the parameters as *{"@S": "name", ...}* with
1147  the remaining list of parameters defined by the *P@S* product model.
[990d8df]1148
1149.. _Test_Your_New_Model:
1150
1151Test Your New Model
1152^^^^^^^^^^^^^^^^^^^
1153
1154Minimal Testing
1155...............
1156
1157From SasView either open the Python shell (*Tools* > *Python Shell/Editor*)
1158or the plugin editor (*Fitting* > *Plugin Model Operations* > *Advanced
1159Plugin Editor*), load your model, and then select *Run > Check Model* from
1160the menu bar. An *Info* box will appear with the results of the compilation
1161and a check that the model runs.
1162
1163Recommended Testing
1164...................
1165
[db1d9d5]1166**NB: For now, this more detailed testing is only possible if you have a
1167SasView build environment available!**
1168
[990d8df]1169If the model compiles and runs, you can next run the unit tests that
1170you have added using the **test =** values.
1171
1172From SasView, switch to the *Shell* tab and type the following::
1173
1174    from sasmodels.model_test import run_one
1175    run_one("~/.sasview/plugin_models/model.py")
1176
1177This should print::
1178
1179    test_model_python (sasmodels.model_test.ModelTestCase) ... ok
1180
1181To check whether single precision is good enough, type the following::
1182
1183    from sasmodels.compare import main as compare
1184    compare("~/.sasview/plugin_models/model.py")
1185
1186This will pop up a plot showing the difference between single precision
1187and double precision on a range of $q$ values.
1188
1189::
1190
1191  demo = dict(scale=1, background=0,
1192              sld=6, sld_solvent=1,
1193              radius=120,
1194              radius_pd=.2, radius_pd_n=45)
1195
1196**demo={'par': value, ...}** in the model file sets the default values for
1197the comparison. You can include polydispersity parameters such as
1198*radius_pd=0.2, radius_pd_n=45* which would otherwise be zero.
1199
1200These commands can also be run directly in the python interpreter:
1201
1202    $ python -m sasmodels.model_test -v ~/.sasview/plugin_models/model.py
1203    $ python -m sasmodels.compare ~/.sasview/plugin_models/model.py
1204
1205The options to compare are quite extensive; type the following for help::
1206
1207    compare()
1208
1209Options will need to be passed as separate strings.
1210For example to run your model with a random set of parameters::
1211
1212    compare("-random", "-pars", "~/.sasview/plugin_models/model.py")
1213
1214For the random models,
1215
1216- *sld* will be in the range (-0.5,10.5),
1217- angles (*theta, phi, psi*) will be in the range (-180,180),
1218- angular dispersion will be in the range (0,45),
1219- polydispersity will be in the range (0,1)
1220- other values will be in the range (0, 2\ *v*), where *v* is the value
1221  of the parameter in demo.
1222
1223Dispersion parameters *n*\, *sigma* and *type* will be unchanged from
1224demo so that run times are more predictable (polydispersity calculated
1225across multiple parameters can be very slow).
1226
[3048ec6]1227If your model has 2D orientation calculation, then you should also
[990d8df]1228test with::
1229
1230    compare("-2d", "~/.sasview/plugin_models/model.py")
1231
1232Check The Docs
1233^^^^^^^^^^^^^^
1234
1235You can get a rough idea of how the documentation will look using the
1236following::
1237
1238    compare("-help", "~/.sasview/plugin_models/model.py")
1239
1240This does not use the same styling as the rest of the docs, but it will
1241allow you to check that your ReStructuredText and LaTeX formatting.
1242Here are some tools to help with the inevitable syntax errors:
1243
1244- `Sphinx cheat sheet <http://matplotlib.org/sampledoc/cheatsheet.html>`_
1245- `Sphinx Documentation <http://www.sphinx-doc.org/en/stable/>`_
1246- `MathJax <http://www.mathjax.org/>`_
1247- `amsmath <http://www.ams.org/publications/authors/tex/amslatex>`_
1248
1249There is also a neat online WYSIWYG ReStructuredText editor at
1250http://rst.ninjs.org\ .
1251
1252
1253Clean Lint - (Developer Version Only)
1254^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1255
1256**NB: For now we are not providing pylint with the installer version
1257of SasView; so unless you have a SasView build environment available,
1258you can ignore this section!**
1259
1260Run the lint check with::
1261
1262    python -m pylint --rcfile=extra/pylint.rc ~/.sasview/plugin_models/model.py
1263
1264We are not aiming for zero lint just yet, only keeping it to a minimum.
1265For now, don't worry too much about *invalid-name*. If you really want a
1266variable name *Rg* for example because $R_g$ is the right name for the model
1267parameter then ignore the lint errors.  Also, ignore *missing-docstring*
[108e70e]1268for standard model functions *Iq*, *Iqac*, etc.
[990d8df]1269
1270We will have delinting sessions at the SasView Code Camps, where we can
1271decide on standards for model files, parameter names, etc.
1272
1273For now, you can tell pylint to ignore things.  For example, to align your
1274parameters in blocks::
1275
1276    # pylint: disable=bad-whitespace,line-too-long
1277    #   ["name",                  "units", default, [lower, upper], "type", "description"],
1278    parameters = [
1279        ["contrast_factor",       "barns",    10.0,  [-inf, inf], "", "Contrast factor of the polymer"],
1280        ["bjerrum_length",        "Ang",       7.1,  [0, inf],    "", "Bjerrum length"],
1281        ["virial_param",          "1/Ang^2",  12.0,  [-inf, inf], "", "Virial parameter"],
1282        ["monomer_length",        "Ang",      10.0,  [0, inf],    "", "Monomer length"],
1283        ["salt_concentration",    "mol/L",     0.0,  [-inf, inf], "", "Concentration of monovalent salt"],
1284        ["ionization_degree",     "",          0.05, [0, inf],    "", "Degree of ionization"],
1285        ["polymer_concentration", "mol/L",     0.7,  [0, inf],    "", "Polymer molar concentration"],
1286        ]
1287    # pylint: enable=bad-whitespace,line-too-long
1288
1289Don't put in too many pylint statements, though, since they make the code ugly.
1290
1291Share Your Model!
1292^^^^^^^^^^^^^^^^^
1293
1294Once compare and the unit test(s) pass properly and everything is done,
1295consider adding your model to the
1296`Model Marketplace <http://marketplace.sasview.org/>`_ so that others may use it!
1297
1298.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
1299
1300*Document History*
1301
1302| 2016-10-25 Steve King
[c654160]1303| 2017-05-07 Paul Kienzle - Moved from sasview to sasmodels docs
[e5cb3df]1304| 2019-03-28 Paul Kienzle - Update docs for radius_effective and shell_volume
Note: See TracBrowser for help on using the repository browser.