source: sasmodels/doc/guide/plugin.rst @ 30b60d2

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 30b60d2 was 30b60d2, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

allow build of both latex and html

  • Property mode set to 100644
File size: 39.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
80- a **short discription**:
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.
235If it is not provided, it will use the name of the model file, with '_'
236replaced by spaces and the parts capitalized.  So *adsorbed_layer.py* will
237become *Adsorbed Layer*.  The predefined models all use the name of the
238model file as the name of the model, so the default may be changed.
239
240**title = "short description"** is short description of the model which
241is included after the model name in the automatically generated documentation.
242The title can also be used for a tooltip.
243
244**description = """doc string"""** is a longer description of the model. It
245shows up when you press the "Description" button of the SasView FitPage.
246It should give a brief description of the equation and the parameters
247without the need to read the entire model documentation. The triple quotes
248allow you to write the description over multiple lines. Keep the lines
249short since the GUI will wrap each one separately if they are too long.
250**Make sure the parameter names in the description match the model definition!**
251
252**category = "shape:sphere"** defines where the model will appear in the
253model documentation.  In this example, the model will appear alphabetically
254in the list of spheroid models in the *Shape* category.
255
256**single = True** indicates that the model can be run using single
257precision floating point values.  Set it to False if the numerical
258calculation for the model is unstable, which is the case for about 20 of
259the built in models.  It is worthwhile modifying the calculation to support
260single precision, allowing models to run up to 10 times faster.  The
261section `Test_Your_New_Model`_  describes how to compare model values for
262single vs. double precision so you can decide if you need to set
263single to False.
264
265**opencl = False** indicates that the model should not be run using OpenCL.
266This may be because the model definition includes code that cannot be
267compiled for the GPU (for example, goto statements).  It can also be used
268for large models which can't run on most GPUs.  This flag has not been
269used on any of the built in models; models which were failing were
270streamlined so this flag was not necessary.
271
272**structure_factor = True** indicates that the model can be used as a
273structure factor to account for interactions between particles.  See
274`Form_Factors`_ for more details.
275
276Model Parameters
277................
278
279Next comes the parameter table.  For example::
280
281    # pylint: disable=bad-whitespace, line-too-long
282    #   ["name",        "units", default, [min, max], "type",    "description"],
283    parameters = [
284        ["sld",         "1e-6/Ang^2",  1, [-inf, inf], "sld",    "Layer scattering length density"],
285        ["sld_solvent", "1e-6/Ang^2",  6, [-inf, inf], "sld",    "Solvent scattering length density"],
286        ["radius",      "Ang",        50, [0, inf],    "volume", "Sphere radius"],
287    ]
288    # pylint: enable=bad-whitespace, line-too-long
289
290**parameters = [["name", "units", default, [min,max], "type", "tooltip"],...]**
291defines the parameters that form the model.
292
293**Note: The order of the parameters in the definition will be the order of the
294parameters in the user interface and the order of the parameters in Iq(),
295Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are
296implicit to all models, so they do not need to be included in the parameter table.**
297
298- **"name"** is the name of the parameter shown on the FitPage.
299
300  - parameter names should follow the mathematical convention; e.g.,
301    *radius_core* not *core_radius*, or *sld_solvent* not *solvent_sld*.
302
303  - model parameter names should be consistent between different models,
304    so *sld_solvent*, for example, should have exactly the same name
305    in every model.
306
307  - to see all the parameter names currently in use, type the following in the
308    python shell/editor under the Tools menu::
309
310       import sasmodels.list_pars
311       sasmodels.list_pars.list_pars()
312
313    *re-use* as many as possible!!!
314
315  - use "name[n]" for multiplicity parameters, where *n* is the name of
316    the parameter defining the number of shells/layers/segments, etc.
317
318- **"units"** are displayed along with the parameter name
319
320  - every parameter should have units; use "None" if there are no units.
321
322  - **sld's should be given in units of 1e-6/Ang^2, and not simply
323    1/Ang^2 to be consistent with the builtin models.  Adjust your formulas
324    appropriately.**
325
326  - fancy units markup is available for some units, including::
327
328        Ang, 1/Ang, 1/Ang^2, 1e-6/Ang^2, degrees, 1/cm, Ang/cm, g/cm^3, mg/m^2
329
330  - the list of units is defined in the variable *RST_UNITS* within
331    `sasmodels/generate.py <https://github.com/SasView/sasmodels/tree/master/sasmodels/generate.py>`_
332
333    - new units can be added using the macros defined in *doc/rst_prolog*
334      in the sasmodels source.
335    - units should be properly formatted using sub-/super-scripts
336      and using negative exponents instead of the / operator, though
337      the unit name should use the / operator for consistency.
338    - please post a message to the SasView developers mailing list with your changes.
339
340- **default** is the initial value for the parameter.
341
342  - **the parameter default values are used to auto-generate a plot of
343    the model function in the documentation.**
344
345- **[min, max]** are the lower and upper limits on the parameter.
346
347  - lower and upper limits can be any number, or *-inf* or *inf*.
348
349  - the limits will show up as the default limits for the fit making it easy,
350    for example, to force the radius to always be greater than zero.
351
352  - these are hard limits defining the valid range of parameter values;
353    polydisperity distributions will be truncated at the limits.
354
355- **"type"** can be one of: "", "sld", "volume", or "orientation".
356
357  - "sld" parameters can have magnetic moments when fitting magnetic models;
358    depending on the spin polarization of the beam and the $q$ value being
359    examined, the effective sld for that material will be used to compute the
360    scattered intensity.
361
362  - "volume" parameters are passed to Iq(), Iqxy(), and form_volume(), and
363    have polydispersity loops generated automatically.
364
365  - "orientation" parameters are only passed to Iqxy(), and have angular
366    dispersion.
367
368
369Model Computation
370.................
371
372Models can be defined as pure python models, or they can be a mixture of
373python and C models.  C models are run on the GPU if it is available,
374otherwise they are compiled and run on the CPU.
375
376Models are defined by the scattering kernel, which takes a set of parameter
377values defining the shape, orientation and material, and returns the
378expected scattering. Polydispersity and angular dispersion are defined
379by the computational infrastructure.  Any parameters defined as "volume"
380parameters are polydisperse, with polydispersity defined in proportion
381to their value.  "orientation" parameters use angular dispersion defined
382in degrees, and are not relative to the current angle.
383
384Based on a weighting function $G(x)$ and a number of points $n$, the
385computed value is
386
387.. math::
388
389     \hat I(q)
390     = \frac{\int G(x) I(q, x)\,dx}{\int G(x) V(x)\,dx}
391     \approx \frac{\sum_{i=1}^n G(x_i) I(q,x_i)}{\sum_{i=1}^n G(x_i) V(x_i)}
392
393That is, the indivdual models do not need to include polydispersity
394calculations, but instead rely on numerical integration to compute the
395appropriately smeared pattern.   Angular dispersion values over polar angle
396$\theta$ requires an additional $\cos \theta$ weighting due to decreased
397arc length for the equatorial angle $\phi$ with increasing latitude.
398
399Python Models
400.............
401
402For pure python models, define the *Iq* function::
403
404      import numpy as np
405      from numpy import cos, sin, ...
406
407      def Iq(q, par1, par2, ...):
408          return I(q, par1, par2, ...)
409      Iq.vectorized = True
410
411The parameters *par1, par2, ...* are the list of non-orientation parameters
412to the model in the order that they appear in the parameter table.
413**Note that the autogenerated model file uses** *x* **rather than** *q*.
414
415The *.py* file should import trigonometric and exponential functions from
416numpy rather than from math.  This lets us evaluate the model for the whole
417range of $q$ values at once rather than looping over each $q$ separately in
418python.  With $q$ as a vector, you cannot use if statements, but must instead
419do tricks like
420
421::
422
423     a = x*q*(q>0) + y*q*(q<=0)
424
425or
426
427::
428
429     a = np.empty_like(q)
430     index = q>0
431     a[index] = x*q[index]
432     a[~index] = y*q[~index]
433
434which sets $a$ to $q \cdot x$ if $q$ is positive or $q \cdot y$ if $q$
435is zero or negative. If you have not converted your function to use $q$
436vectors, you can set the following and it will only receive one $q$
437value at a time::
438
439    Iq.vectorized = False
440
441Return np.NaN if the parameters are not valid (e.g., cap_radius < radius in
442barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be
443ignored, and not included in the calculation of the weighted polydispersity.
444
445Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the
446parameter list includes any orientation parameters.  If *Iqxy* is not defined,
447then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*.
448
449Models should define *form_volume(par1, par2, ...)* where the parameter
450list includes the *volume* parameters in order.  This is used for a weighted
451volume normalization so that scattering is on an absolute scale.  If
452*form_volume* is not defined, then the default *form_volume = 1.0* will be
453used.
454
455Embedded C Models
456.................
457
458Like pure python models, inline C models need to define an *Iq* function::
459
460    Iq = """
461        return I(q, par1, par2, ...);
462    """
463
464This expands into the equivalent C code::
465
466    #include <math.h>
467    double Iq(double q, double par1, double par2, ...);
468    double Iq(double q, double par1, double par2, ...)
469    {
470        return I(q, par1, par2, ...);
471    }
472
473*Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*,
474and it includes orientation parameters.
475
476*form_volume* defines the volume of the shape. As in python models, it
477includes only the volume parameters.
478
479*Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and
480*form_volume* will default to 1.0.
481
482**source=['fn.c', ...]** includes the listed C source files in the
483program before *Iq* and *Iqxy* are defined. This allows you to extend the
484library of C functions available to your model.
485
486Models are defined using double precision declarations for the
487parameters and return values.  When a model is run using single
488precision or long double precision, each variable is converted
489to the target type, depending on the precision requested.
490
491**Floating point constants must include the decimal point.**  This allows us
492to convert values such as 1.0 (double precision) to 1.0f (single precision)
493so that expressions that use these values are not promoted to double precision
494expressions.  Some graphics card drivers are confused when functions
495that expect floating point values are passed integers, such as 4*atan(1); it
496is safest to not use integers in floating point expressions.  Even better,
497use the builtin constant M_PI rather than 4*atan(1); it is faster and smaller!
498
499The C model operates on a single $q$ value at a time.  The code will be
500run in parallel across different $q$ values, either on the graphics card
501or the processor.
502
503Rather than returning NAN from Iq, you must define the *INVALID(v)*.  The
504*v* parameter lets you access all the parameters in the model using
505*v.par1*, *v.par2*, etc. For example::
506
507    #define INVALID(v) (v.bell_radius < v.radius)
508
509Special Functions
510.................
511
512The C code follows the C99 standard, with the usual math functions,
513as defined in
514`OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_.
515This includes the following:
516
517    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E:
518        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$
519    exp, log, pow(x,y), expm1, sqrt:
520        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$.
521        The function expm1(x) is accurate across all $x$, including $x$
522        very close to zero.
523    sin, cos, tan, asin, acos, atan:
524        Trigonometry functions and inverses, operating on radians.
525    sinh, cosh, tanh, asinh, acosh, atanh:
526        Hyperbolic trigonometry functions.
527    atan2(y,x):
528        Angle from the $x$\ -axis to the point $(x,y)$, which is equal to
529        $\tan^{-1}(y/x)$ corrected for quadrant.  That is, if $x$ and $y$ are
530        both negative, then atan2(y,x) returns a value in quadrant III where
531        atan(y/x) would return a value in quadrant I. Similarly for
532        quadrants II and IV when $x$ and $y$ have opposite sign.
533    fmin(x,y), fmax(x,y), trunc, rint:
534        Floating point functions.  rint(x) returns the nearest integer.
535    NAN:
536        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that
537        you cannot use :code:`x == NAN` to test for NaN values since that
538        will always return false.  NAN does not equal NAN!
539    INFINITY:
540        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x)
541        to test for finite and not NaN.
542    erf, erfc, tgamma, lgamma:  **do not use**
543        Special functions that should be part of the standard, but are missing
544        or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma
545        instead (see below). Note: lgamma(x) has not yet been tested.
546
547Some non-standard constants and functions are also provided:
548
549    M_PI_180, M_4PI_3:
550        $\frac{\pi}{180}$, $\frac{4\pi}{3}$
551    SINCOS(x, s, c):
552        Macro which sets s=sin(x) and c=cos(x). The variables *c* and *s*
553        must be declared first.
554    square(x):
555        $x^2$
556    cube(x):
557        $x^3$
558    sas_sinx_x(x):
559        $\sin(x)/x$, with limit $\sin(0)/0 = 1$.
560    powr(x, y):
561        $x^y$ for $x \ge 0$; this is faster than general $x^y$ on some GPUs.
562    pown(x, n):
563        $x^n$ for $n$ integer; this is faster than general $x^n$ on some GPUs.
564    FLOAT_SIZE:
565        The number of bytes in a floating point value.  Even though all
566        variables are declared double, they may be converted to single
567        precision float before running. If your algorithm depends on
568        precision (which is not uncommon for numerical algorithms), use
569        the following::
570
571            #if FLOAT_SIZE>4
572            ... code for double precision ...
573            #else
574            ... code for single precision ...
575            #endif
576    SAS_DOUBLE:
577        A replacement for :code:`double` so that the declared variable will
578        stay double precision; this should generally not be used since some
579        graphics cards do not support double precision.  There is no provision
580        for forcing a constant to stay double precision.
581
582The following special functions and scattering calculations are defined in
583`sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_.
584These functions have been tuned to be fast and numerically stable down
585to $q=0$ even in single precision.  In some cases they work around bugs
586which appear on some platforms but not others, so use them where needed.
587Add the files listed in :code:`source = ["lib/file.c", ...]` to your *model.py*
588file in the order given, otherwise these functions will not be available.
589
590    polevl(x, c, n):
591        Polynomial evaluation $p(x) = \sum_{i=0}^n c_i x^i$ using Horner's
592        method so it is faster and more accurate.
593
594        $c = \{c_n, c_{n-1}, \ldots, c_0 \}$ is the table of coefficients,
595        sorted from highest to lowest.
596
597        :code:`source = ["lib/polevl.c", ...]` (`link to code <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
598
599    p1evl(x, c, n):
600        Evaluation of normalized polynomial $p(x) = x^n + \sum_{i=0}^{n-1} c_i x^i$
601        using Horner's method so it is faster and more accurate.
602
603        $c = \{c_{n-1}, c_{n-2} \ldots, c_0 \}$ is the table of coefficients,
604        sorted from highest to lowest.
605
606        :code:`source = ["lib/polevl.c", ...]`
[870a2f4]607        (`polevl.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/polevl.c>`_)
[990d8df]608
609    sas_gamma(x):
[30b60d2]610        Gamma function sas_gamma\ $(x) = \Gamma(x)$.
[990d8df]611
612        The standard math function, tgamma(x) is unstable for $x < 1$
613        on some platforms.
614
[870a2f4]615        :code:`source = ["lib/sas_gamma.c", ...]`
616        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_)
[990d8df]617
618    sas_erf(x), sas_erfc(x):
619        Error function
[30b60d2]620        sas_erf\ $(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt$
[990d8df]621        and complementary error function
[30b60d2]622        sas_erfc\ $(x) = \frac{2}{\sqrt\pi}\int_x^{\infty} e^{-t^2}\,dt$.
[990d8df]623
624        The standard math functions erf(x) and erfc(x) are slower and broken
625        on some platforms.
626
627        :code:`source = ["lib/polevl.c", "lib/sas_erf.c", ...]`
[870a2f4]628        (`sas_erf.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_erf.c>`_)
[990d8df]629
630    sas_J0(x):
[30b60d2]631        Bessel function of the first kind sas_J0\ $(x)=J_0(x)$ where
[990d8df]632        $J_0(x) = \frac{1}{\pi}\int_0^\pi \cos(x\sin(\tau))\,d\tau$.
633
634        The standard math function j0(x) is not available on all platforms.
635
636        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", ...]`
[870a2f4]637        (`sas_J0.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J0.c>`_)
[990d8df]638
639    sas_J1(x):
[30b60d2]640        Bessel function of the first kind  sas_J1\ $(x)=J_1(x)$ where
[990d8df]641        $J_1(x) = \frac{1}{\pi}\int_0^\pi \cos(\tau - x\sin(\tau))\,d\tau$.
642
643        The standard math function j1(x) is not available on all platforms.
644
645        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
[870a2f4]646        (`sas_J1.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
[990d8df]647
648    sas_JN(n, x):
[30b60d2]649        Bessel function of the first kind and integer order $n$,
650        sas_JN\ $(n, x) =J_n(x)$ where
[990d8df]651        $J_n(x) = \frac{1}{\pi}\int_0^\pi \cos(n\tau - x\sin(\tau))\,d\tau$.
[30b60d2]652        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively.
[990d8df]653
654        The standard math function jn(n, x) is not available on all platforms.
655
656        :code:`source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", "lib/sas_JN.c", ...]`
[870a2f4]657        (`sas_JN.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_JN.c>`_)
[990d8df]658
659    sas_Si(x):
[30b60d2]660        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
[990d8df]661
662        This function uses Taylor series for small and large arguments:
663
664        For large arguments,
665
666        .. math::
667
668             \text{Si}(x) \sim \frac{\pi}{2}
669             - \frac{\cos(x)}{x}\left(1 - \frac{2!}{x^2} + \frac{4!}{x^4} - \frac{6!}{x^6} \right)
670             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)
671
672        For small arguments,
673
674        .. math::
675
676           \text{Si}(x) \sim x
677           - \frac{x^3}{3\times 3!} + \frac{x^5}{5 \times 5!} - \frac{x^7}{7 \times 7!}
678           + \frac{x^9}{9\times 9!} - \frac{x^{11}}{11\times 11!}
679
680        :code:`source = ["lib/Si.c", ...]`
[870a2f4]681        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_)
[990d8df]682
683    sas_3j1x_x(x):
684        Spherical Bessel form
[30b60d2]685        sph_j1c\ $(x) = 3 j_1(x)/x = 3 (\sin(x) - x \cos(x))/x^3$,
[990d8df]686        with a limiting value of 1 at $x=0$, where $j_1(x)$ is the spherical
687        Bessel function of the first kind and first order.
688
689        This function uses a Taylor series for small $x$ for numerical accuracy.
690
691        :code:`source = ["lib/sas_3j1x_x.c", ...]`
[870a2f4]692        (`sas_3j1x_x.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_3j1x_x.c>`_)
[990d8df]693
694
695    sas_2J1x_x(x):
[30b60d2]696        Bessel form sas_J1c\ $(x) = 2 J_1(x)/x$, with a limiting value
[990d8df]697        of 1 at $x=0$, where $J_1(x)$ is the Bessel function of first kind
698        and first order.
699
700        :code:`source = ["lib/polevl.c", "lib/sas_J1.c", ...]`
[870a2f4]701        (`sas_J1.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_J1.c>`_)
[990d8df]702
703
704    Gauss76Z[i], Gauss76Wt[i]:
705        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively,
706        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$.
707
708        Similar arrays are available in :code:`gauss20.c` for 20-point
709        quadrature and in :code:`gauss150.c` for 150-point quadrature.
710
711        :code:`source = ["lib/gauss76.c", ...]`
[870a2f4]712        (`gauss76.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/gauss76.c>`_)
[990d8df]713
714
715
716Problems with C models
717......................
718
719The graphics processor (GPU) in your computer is a specialized computer tuned
720for certain kinds of problems.  This leads to strange restrictions that you
721need to be aware of.  Your code may work fine on some platforms or for some
722models, but then return bad values on other platforms.  Some examples of
723particular problems:
724
725  **(1) Code is too complex, or uses too much memory.** GPU devices only
726  have a limited amount of memory available for each processor. If you run
727  programs which take too much memory, then rather than running multiple
728  values in parallel as it usually does, the GPU may only run a single
729  version of the code at a time, making it slower than running on the CPU.
730  It may fail to run on some platforms, or worse, cause the screen to go
731  blank or the system to reboot.
732
733  **(2) Code takes too long.** Because GPU devices are used for the computer
734  display, the OpenCL drivers are very careful about the amount of time they
735  will allow any code to run. For example, on OS X, the model will stop
736  running after 5 seconds regardless of whether the computation is complete.
737  You may end up with only some of your 2D array defined, with the rest
738  containing random data. Or it may cause the screen to go blank or the
739  system to reboot.
740
741  **(3) Memory is not aligned**. The GPU hardware is specialized to operate
742  on multiple values simultaneously. To keep the GPU simple the values in
743  memory must be aligned with the different GPU compute engines. Not
744  following these rules can lead to unexpected values being loaded into
745  memory, and wrong answers computed. The conclusion from a very long and
746  strange debugging session was that any arrays that you declare in your
747  model should be a multiple of four. For example::
748
749      double Iq(q, p1, p2, ...)
750      {
751          double vector[8];  // Only going to use seven slots, but declare 8
752          ...
753      }
754
755The first step when your model is behaving strangely is to set
756**single=False**. This automatically restricts the model to only run on the
757CPU, or on high-end GPU cards. There can still be problems even on high-end
758cards, so you can force the model off the GPU by setting **opencl=False**.
759This runs the model as a normal C program without any GPU restrictions so
760you know that strange results are probably from your code rather than the
761environment. Once the code is debugged, you can compare your output to the
762output on the GPU.
763
764Although it can be difficult to get your model to work on the GPU, the reward
765can be a model that runs 1000x faster on a good card.  Even your laptop may
766show a 50x improvement or more over the equivalent pure python model.
767
768External C Models
769.................
770
771External C models are very much like embedded C models, except that
772*Iq*, *Iqxy* and *form_volume* are defined in an external source file
773loaded using the *source=[...]* statement. You need to supply the function
774declarations for each of these that you need instead of building them
775automatically from the parameter table.
776
777
778.. _Form_Factors:
779
780Form Factors
781............
782
783Away from the dilute limit you can estimate scattering including
784particle-particle interactions using $I(q) = P(q)*S(q)$ where $P(q)$
785is the form factor and $S(q)$ is the structure factor.  The simplest
786structure factor is the *hardsphere* interaction, which
787uses the effective radius of the form factor as an input to the structure
788factor model.  The effective radius is the average radius of the
789form averaged over all the polydispersity values.
790
791::
792
793    def ER(radius, thickness):
794        """Effective radius of a core-shell sphere."""
795        return radius + thickness
796
797Now consider the *core_shell_sphere*, which has a simple effective radius
798equal to the radius of the core plus the thickness of the shell, as
799shown above. Given polydispersity over *(r1, r2, ..., rm)* in radius and
800*(t1, t2, ..., tn)* in thickness, *ER* is called with a mesh
801grid covering all possible combinations of radius and thickness.
802That is, *radius* is *(r1, r2, ..., rm, r1, r2, ..., rm, ...)*
803and *thickness* is *(t1, t1, ... t1, t2, t2, ..., t2, ...)*.
804The *ER* function returns one effective radius for each combination.
805The effective radius calculator weights each of these according to
806the polydispersity distributions and calls the structure factor
807with the average *ER*.
808
809::
810
811    def VR(radius, thickness):
812        """Sphere and shell volumes for a core-shell sphere."""
813        whole = 4.0/3.0 * pi * (radius + thickness)**3
814        core = 4.0/3.0 * pi * radius**3
815        return whole, whole - core
816
817Core-shell type models have an additional volume ratio which scales
818the structure factor.  The *VR* function returns the volume of
819the whole sphere and the volume of the shell. Like *ER*, there is
820one return value for each point in the mesh grid.
821
822*NOTE: we may be removing or modifying this feature soon. As of the
823time of writing, core-shell sphere returns (1., 1.) for VR, giving a volume
824ratio of 1.0.*
825
826Unit Tests
827..........
828
829THESE ARE VERY IMPORTANT. Include at least one test for each model and
830PLEASE make sure that the answer value is correct (i.e. not a random number).
831
832::
833
834    tests = [
835        [{}, 0.2, 0.726362],
836        [{"scale": 1., "background": 0., "sld": 6., "sld_solvent": 1.,
837          "radius": 120., "radius_pd": 0.2, "radius_pd_n":45},
838         0.2, 0.228843],
839        [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.],
840        [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.],
841    ]
842
843
844**tests=[[{parameters}, q, result], ...]** is a list of lists.
845Each list is one test and contains, in order:
846
847- a dictionary of parameter values. This can be *{}* using the default
848  parameters, or filled with some parameters that will be different from the
849  default, such as *{"radius":10.0, "sld":4}*. Unlisted parameters will
850  be given the default values.
851- the input $q$ value or tuple of $(q_x, q_y)$ values.
852- the output $I(q)$ or $I(q_x,q_y)$ expected of the model for the parameters
853  and input value given.
854- input and output values can themselves be lists if you have several
855  $q$ values to test for the same model parameters.
856- for testing *ER* and *VR*, give the inputs as "ER" and "VR" respectively;
857  the output for *VR* should be the sphere/shell ratio, not the individual
858  sphere and shell values.
859
860.. _Test_Your_New_Model:
861
862Test Your New Model
863^^^^^^^^^^^^^^^^^^^
864
865Minimal Testing
866...............
867
868From SasView either open the Python shell (*Tools* > *Python Shell/Editor*)
869or the plugin editor (*Fitting* > *Plugin Model Operations* > *Advanced
870Plugin Editor*), load your model, and then select *Run > Check Model* from
871the menu bar. An *Info* box will appear with the results of the compilation
872and a check that the model runs.
873
874If you are not using sasmodels from SasView, skip this step.
875
876Recommended Testing
877...................
878
879If the model compiles and runs, you can next run the unit tests that
880you have added using the **test =** values.
881
882From SasView, switch to the *Shell* tab and type the following::
883
884    from sasmodels.model_test import run_one
885    run_one("~/.sasview/plugin_models/model.py")
886
887This should print::
888
889    test_model_python (sasmodels.model_test.ModelTestCase) ... ok
890
891To check whether single precision is good enough, type the following::
892
893    from sasmodels.compare import main as compare
894    compare("~/.sasview/plugin_models/model.py")
895
896This will pop up a plot showing the difference between single precision
897and double precision on a range of $q$ values.
898
899::
900
901  demo = dict(scale=1, background=0,
902              sld=6, sld_solvent=1,
903              radius=120,
904              radius_pd=.2, radius_pd_n=45)
905
906**demo={'par': value, ...}** in the model file sets the default values for
907the comparison. You can include polydispersity parameters such as
908*radius_pd=0.2, radius_pd_n=45* which would otherwise be zero.
909
910These commands can also be run directly in the python interpreter:
911
912    $ python -m sasmodels.model_test -v ~/.sasview/plugin_models/model.py
913    $ python -m sasmodels.compare ~/.sasview/plugin_models/model.py
914
915The options to compare are quite extensive; type the following for help::
916
917    compare()
918
919Options will need to be passed as separate strings.
920For example to run your model with a random set of parameters::
921
922    compare("-random", "-pars", "~/.sasview/plugin_models/model.py")
923
924For the random models,
925
926- *sld* will be in the range (-0.5,10.5),
927- angles (*theta, phi, psi*) will be in the range (-180,180),
928- angular dispersion will be in the range (0,45),
929- polydispersity will be in the range (0,1)
930- other values will be in the range (0, 2\ *v*), where *v* is the value
931  of the parameter in demo.
932
933Dispersion parameters *n*\, *sigma* and *type* will be unchanged from
934demo so that run times are more predictable (polydispersity calculated
935across multiple parameters can be very slow).
936
937If your model has 2D orientational calculation, then you should also
938test with::
939
940    compare("-2d", "~/.sasview/plugin_models/model.py")
941
942Check The Docs
943^^^^^^^^^^^^^^
944
945You can get a rough idea of how the documentation will look using the
946following::
947
948    compare("-help", "~/.sasview/plugin_models/model.py")
949
950This does not use the same styling as the rest of the docs, but it will
951allow you to check that your ReStructuredText and LaTeX formatting.
952Here are some tools to help with the inevitable syntax errors:
953
954- `Sphinx cheat sheet <http://matplotlib.org/sampledoc/cheatsheet.html>`_
955- `Sphinx Documentation <http://www.sphinx-doc.org/en/stable/>`_
956- `MathJax <http://www.mathjax.org/>`_
957- `amsmath <http://www.ams.org/publications/authors/tex/amslatex>`_
958
959There is also a neat online WYSIWYG ReStructuredText editor at
960http://rst.ninjs.org\ .
961
962
963Clean Lint - (Developer Version Only)
964^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
965
966**NB: For now we are not providing pylint with the installer version
967of SasView; so unless you have a SasView build environment available,
968you can ignore this section!**
969
970Run the lint check with::
971
972    python -m pylint --rcfile=extra/pylint.rc ~/.sasview/plugin_models/model.py
973
974We are not aiming for zero lint just yet, only keeping it to a minimum.
975For now, don't worry too much about *invalid-name*. If you really want a
976variable name *Rg* for example because $R_g$ is the right name for the model
977parameter then ignore the lint errors.  Also, ignore *missing-docstring*
978for standard model functions *Iq*, *Iqxy*, etc.
979
980We will have delinting sessions at the SasView Code Camps, where we can
981decide on standards for model files, parameter names, etc.
982
983For now, you can tell pylint to ignore things.  For example, to align your
984parameters in blocks::
985
986    # pylint: disable=bad-whitespace,line-too-long
987    #   ["name",                  "units", default, [lower, upper], "type", "description"],
988    parameters = [
989        ["contrast_factor",       "barns",    10.0,  [-inf, inf], "", "Contrast factor of the polymer"],
990        ["bjerrum_length",        "Ang",       7.1,  [0, inf],    "", "Bjerrum length"],
991        ["virial_param",          "1/Ang^2",  12.0,  [-inf, inf], "", "Virial parameter"],
992        ["monomer_length",        "Ang",      10.0,  [0, inf],    "", "Monomer length"],
993        ["salt_concentration",    "mol/L",     0.0,  [-inf, inf], "", "Concentration of monovalent salt"],
994        ["ionization_degree",     "",          0.05, [0, inf],    "", "Degree of ionization"],
995        ["polymer_concentration", "mol/L",     0.7,  [0, inf],    "", "Polymer molar concentration"],
996        ]
997    # pylint: enable=bad-whitespace,line-too-long
998
999Don't put in too many pylint statements, though, since they make the code ugly.
1000
1001Share Your Model!
1002^^^^^^^^^^^^^^^^^
1003
1004Once compare and the unit test(s) pass properly and everything is done,
1005consider adding your model to the
1006`Model Marketplace <http://marketplace.sasview.org/>`_ so that others may use it!
1007
1008.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
1009
1010*Document History*
1011
1012| 2016-10-25 Steve King
1013| 2017-05-07 Paul Kienzle - Moved from sasview to sasmodels docs
Note: See TracBrowser for help on using the repository browser.