source: sasmodels/sasmodels/generate.py @ 4eac427

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 4eac427 was 4eac427, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

fix doc build for units without special markup

  • Property mode set to 100644
File size: 22.6 KB
Line 
1"""
2SAS model constructor.
3
4Small angle scattering models are defined by a set of kernel functions:
5
6    *Iq(q, p1, p2, ...)* returns the scattering at q for a form with
7    particular dimensions averaged over all orientations.
8
9    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx,qy for a form
10    with particular dimensions for a single orientation.
11
12    *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the
13    polarized neutron spin states (up-up, up-down, down-up, down-down) for
14    a form with particular dimensions for a single orientation.
15
16    *form_volume(p1, p2, ...)* returns the volume of the form with particular
17    dimension.
18
19    *ER(p1, p2, ...)* returns the effective radius of the form with
20    particular dimensions.
21
22    *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms.
23
24These functions are defined in a kernel module .py script and an associated
25set of .c files.  The model constructor will use them to create models with
26polydispersity across volume and orientation parameters, and provide
27scale and background parameters for each model.
28
29*Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99
30functions written for OpenCL.  All functions need prototype declarations
31even if the are defined before they are used.  OpenCL does not support
32*#include* preprocessor directives, so instead the list of includes needs
33to be given as part of the metadata in the kernel module definition.
34The included files should be listed using a path relative to the kernel
35module, or if using "lib/file.c" if it is one of the standard includes
36provided with the sasmodels source.  The includes need to be listed in
37order so that functions are defined before they are used.
38
39Floating point values should be declared as *double*.  For single precision
40calculations, *double* will be replaced by *float*.  The single precision
41conversion will also tag floating point constants with "f" to make them
42single precision constants.  When using integral values in floating point
43expressions, they should be expressed as floating point values by including
44a decimal point.  This includes 0., 1. and 2.
45
46OpenCL has a *sincos* function which can improve performance when both
47the *sin* and *cos* values are needed for a particular argument.  Since
48this function does not exist in C99, all use of *sincos* should be
49replaced by the macro *SINCOS(value,sn,cn)* where *sn* and *cn* are
50previously declared *double* variables.  When compiled for systems without
51OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.   If *value* is
52an expression, it will appear twice in this case; whether or not it will be
53evaluated twice depends on the quality of the compiler.
54
55If the input parameters are invalid, the scattering calculator should
56return a negative number. Particularly with polydispersity, there are
57some sets of shape parameters which lead to nonsensical forms, such
58as a capped cylinder where the cap radius is smaller than the
59cylinder radius.  The polydispersity calculation will ignore these points,
60effectively chopping the parameter weight distributions at the boundary
61of the infeasible region.  The resulting scattering will be set to
62background.  This will work correctly even when polydispersity is off.
63
64*ER* and *VR* are python functions which operate on parameter vectors.
65The constructor code will generate the necessary vectors for computing
66them with the desired polydispersity.
67
68The available kernel parameters are defined as a list, with each parameter
69defined as a sublist with the following elements:
70
71    *name* is the name that will be used in the call to the kernel
72    function and the name that will be displayed to the user.  Names
73    should be lower case, with words separated by underscore.  If
74    acronyms are used, the whole acronym should be upper case.
75
76    *units* should be one of *degrees* for angles, *Ang* for lengths,
77    *1e-6/Ang^2* for SLDs.
78
79    *default value* will be the initial value for  the model when it
80    is selected, or when an initial value is not otherwise specified.
81
82    [*lb*, *ub*] are the hard limits on the parameter value, used to limit
83    the polydispersity density function.  In the fit, the parameter limits
84    given to the fit are the limits  on the central value of the parameter.
85    If there is polydispersity, it will evaluate parameter values outside
86    the fit limits, but not outside the hard limits specified in the model.
87    If there are no limits, use +/-inf imported from numpy.
88
89    *type* indicates how the parameter will be used.  "volume" parameters
90    will be used in all functions.  "orientation" parameters will be used
91    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
92    *Imagnetic* only.  If *type* is the empty string, the parameter will
93    be used in all of *Iq*, *Iqxy* and *Imagnetic*.
94
95    *description* is a short description of the parameter.  This will
96    be displayed in the parameter table and used as a tool tip for the
97    parameter value in the user interface.
98
99The kernel module must set variables defining the kernel meta data:
100
101    *name* is the model name
102
103    *title* is a short description of the model, suitable for a tool tip,
104    or a one line model summary in a table of models.
105
106    *description* is an extended description of the model to be displayed
107    while the model parameters are being edited.
108
109    *parameters* is the list of parameters.  Parameters in the kernel
110    functions must appear in the same order as they appear in the
111    parameters list.  Two additional parameters, *scale* and *background*
112    are added to the beginning of the parameter list.  They will show up
113    in the documentation as model parameters, but they are never sent to
114    the kernel functions.
115
116    *source* is the list of C-99 source files that must be joined to
117    create the OpenCL kernel functions.  The files defining the functions
118    need to be listed before the files which use the functions.
119
120    *ER* is a python function defining the effective radius.  If it is
121    not present, the effective radius is 0.
122
123    *VR* is a python function defining the volume ratio.  If it is not
124    present, the volume ratio is 1.
125
126    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the
127    C source code for the body of the volume, Iq, and Iqxy functions
128    respectively.  These can also be defined in the last source file.
129
130    *Iq* and *Iqxy* also be instead be python functions defining the
131    kernel.  If they are marked as *Iq.vectorized = True* then the
132    kernel is passed the entire *q* vector at once, otherwise it is
133    passed values one *q* at a time.  The performance improvement of
134    this step is significant.
135
136An *info* dictionary is constructed from the kernel meta data and
137returned to the caller.
138
139Additional fields can be defined in the kernel definition file that
140are not needed for sas modelling.
141
142    *demo* is a dictionary of parameter=value defining a set of
143    parameters to use by default when *compare* is called.
144
145    *oldname* is the name of the model in sasview before sasmodels
146    was split into its own package, and *oldpars* is a dictionary
147    of *parameter: old_parameter* pairs defining the new names for
148    the parameters.  This is used by *compare* to check the values
149    of the new model against the values of the old model before
150    you are ready to add the new model to sasmodels.
151
152The model evaluator, function call sequence consists of q inputs and the return vector,
153followed by the loop value/weight vector, followed by the values for
154the non-polydisperse parameters, followed by the lengths of the
155polydispersity loops.  To construct the call for 1D models, the
156categories *fixed-1d* and *pd-1d* list the names of the parameters
157of the non-polydisperse and the polydisperse parameters respectively.
158Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models.
159The *pd-rel* category is a set of those parameters which give
160polydispersitiy as a portion of the value (so a 10% length dispersity
161would use a polydispersity value of 0.1) rather than absolute
162dispersity such as an angle plus or minus 15 degrees.
163
164The *volume* category lists the volume parameters in order for calls
165to volume within the kernel (used for volume normalization) and for
166calls to ER and VR for effective radius and volume ratio respectively.
167
168The *orientation* and *magnetic* categories list the orientation and
169magnetic parameters.  These are used by the sasview interface.  The
170blank category is for parameters such as scale which don't have any
171other marking.
172
173The doc string at the start of the kernel module will be used to
174construct the model documentation web pages.  Embedded figures should
175appear in the subdirectory "img" beside the model definition, and tagged
176with the kernel module name to avoid collision with other models.  Some
177file systems are case-sensitive, so only use lower case characters for
178file names and extensions.
179
180
181The function :func:`make` loads the metadata from the module and returns
182the kernel source.  The function :func:`doc` extracts the doc string
183and adds the parameter table to the top.  The function :func:`sources`
184returns a list of files required by the model.
185"""
186
187# TODO: identify model files which have changed since loading and reload them.
188
189__all__ = ["make", "doc", "sources", "use_single"]
190
191import sys
192from os.path import abspath, dirname, join as joinpath, exists
193import re
194
195import numpy as np
196C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
197
198F64 = np.dtype('float64')
199F32 = np.dtype('float32')
200
201# Scale and background, which are parameters common to every form factor
202COMMON_PARAMETERS = [
203    ["scale", "", 1, [0, np.inf], "", "Source intensity"],
204    ["background", "1/cm", 0, [0, np.inf], "", "Source background"],
205    ]
206
207
208# Conversion from units defined in the parameter table for each model
209# to units displayed in the sphinx documentation.
210RST_UNITS = {
211    "Ang": "|Ang|",
212    "1/Ang": "|Ang^-1|",
213    "1/Ang^2": "|Ang^-2|",
214    "1e-6/Ang^2": "|1e-6Ang^-2|",
215    "degrees": "degree",
216    "1/cm": "|cm^-1|",
217    "": "None",
218    }
219
220# Headers for the parameters tables in th sphinx documentation
221PARTABLE_HEADERS = [
222    "Parameter",
223    "Description",
224    "Units",
225    "Default value",
226    ]
227
228# Minimum width for a default value (this is shorter than the column header
229# width, so will be ignored).
230PARTABLE_VALUE_WIDTH = 10
231
232# Documentation header for the module, giving the model name, its short
233# description and its parameter table.  The remainder of the doc comes
234# from the module docstring.
235DOC_HEADER = """.. _%(name)s:
236
237%(label)s
238=======================================================
239
240%(title)s
241
242%(parameters)s
243
244The returned value is scaled to units of |cm^-1|.
245
246%(docs)s
247"""
248
249def format_units(par):
250    return RST_UNITS.get(par, par)
251
252def make_partable(pars):
253    """
254    Generate the parameter table to include in the sphinx documentation.
255    """
256    pars = COMMON_PARAMETERS + pars
257    column_widths = [
258        max(len(p[0]) for p in pars),
259        max(len(p[-1]) for p in pars),
260        max(len(format_units(p[1])) for p in pars),
261        PARTABLE_VALUE_WIDTH,
262        ]
263    column_widths = [max(w, len(h))
264                     for w, h in zip(column_widths, PARTABLE_HEADERS)]
265
266    sep = " ".join("="*w for w in column_widths)
267    lines = [
268        sep,
269        " ".join("%-*s" % (w, h) for w, h in zip(column_widths, PARTABLE_HEADERS)),
270        sep,
271        ]
272    for p in pars:
273        lines.append(" ".join([
274            "%-*s" % (column_widths[0], p[0]),
275            "%-*s" % (column_widths[1], p[-1]),
276            "%-*s" % (column_widths[2], format_units(p[1])),
277            "%*g" % (column_widths[3], p[2]),
278            ]))
279    lines.append(sep)
280    return "\n".join(lines)
281
282def _search(search_path, filename):
283    """
284    Find *filename* in *search_path*.
285
286    Raises ValueError if file does not exist.
287    """
288    for path in search_path:
289        target = joinpath(path, filename)
290        if exists(target):
291            return target
292    raise ValueError("%r not found in %s" % (filename, search_path))
293
294def sources(info):
295    """
296    Return a list of the sources file paths for the module.
297    """
298    search_path = [dirname(info['filename']),
299                   abspath(joinpath(dirname(__file__), 'models'))]
300    return [_search(search_path, f) for f in info['source']]
301
302def use_single(source):
303    """
304    Convert code from double precision to single precision.
305    """
306    # Convert double keyword to float.  Accept an 'n' parameter for vector
307    # values, where n is 2, 4, 8 or 16. Assume complex numbers are represented
308    # as cdouble which is typedef'd to double2.
309    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
310                    r'\1float\2', source)
311    # Convert floating point constants to single by adding 'f' to the end.
312    # OS/X driver complains if you don't do this.
313    source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?',
314                    r'\g<0>f', source)
315    return source
316
317
318def kernel_name(info, is_2D):
319    """
320    Name of the exported kernel symbol.
321    """
322    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq")
323
324
325def categorize_parameters(pars):
326    """
327    Build parameter categories out of the the parameter definitions.
328
329    Returns a dictionary of categories.
330    """
331    partype = {
332        'volume': [], 'orientation': [], 'magnetic': [], '': [],
333        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [],
334        'pd-rel': set(),
335    }
336
337    for p in pars:
338        name, ptype = p[0], p[4]
339        if ptype == 'volume':
340            partype['pd-1d'].append(name)
341            partype['pd-2d'].append(name)
342            partype['pd-rel'].add(name)
343        elif ptype == 'magnetic':
344            partype['fixed-2d'].append(name)
345        elif ptype == 'orientation':
346            partype['pd-2d'].append(name)
347        elif ptype == '':
348            partype['fixed-1d'].append(name)
349            partype['fixed-2d'].append(name)
350        else:
351            raise ValueError("unknown parameter type %r" % ptype)
352        partype[ptype].append(name)
353
354    return partype
355
356def indent(s, depth):
357    """
358    Indent a string of text with *depth* additional spaces on each line.
359    """
360    spaces = " "*depth
361    sep = "\n" + spaces
362    return spaces + sep.join(s.split("\n"))
363
364
365def build_polydispersity_loops(pd_pars):
366    """
367    Build polydispersity loops
368
369    Returns loop opening and loop closing
370    """
371    LOOP_OPEN = """\
372for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) {
373  const double %(name)s = loops[2*(%(name)s_i%(offset)s)];
374  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\
375"""
376    depth = 4
377    offset = ""
378    loop_head = []
379    loop_end = []
380    for name in pd_pars:
381        subst = {'name': name, 'offset': offset}
382        loop_head.append(indent(LOOP_OPEN % subst, depth))
383        loop_end.insert(0, (" "*depth) + "}")
384        offset += '+N' + name
385        depth += 2
386    return "\n".join(loop_head), "\n".join(loop_end)
387
388C_KERNEL_TEMPLATE = None
389def make_model(info):
390    """
391    Generate the code for the kernel defined by info, using source files
392    found in the given search path.
393    """
394    # TODO: need something other than volume to indicate dispersion parameters
395    # No volume normalization despite having a volume parameter.
396    # Thickness is labelled a volume in order to trigger polydispersity.
397    # May want a separate dispersion flag, or perhaps a separate category for
398    # disperse, but not volume.  Volume parameters also use relative values
399    # for the distribution rather than the absolute values used by angular
400    # dispersion.  Need to be careful that necessary parameters are available
401    # for computing volume even if we allow non-disperse volume parameters.
402
403    # Load template
404    global C_KERNEL_TEMPLATE
405    if C_KERNEL_TEMPLATE is None:
406        with open(C_KERNEL_TEMPLATE_PATH) as fid:
407            C_KERNEL_TEMPLATE = fid.read()
408
409    # Load additional sources
410    source = [open(f).read() for f in sources(info)]
411
412    # Prepare defines
413    defines = []
414    partype = info['partype']
415    pd_1d = partype['pd-1d']
416    pd_2d = partype['pd-2d']
417    fixed_1d = partype['fixed-1d']
418    fixed_2d = partype['fixed-1d']
419
420    iq_parameters = [p[0]
421                     for p in info['parameters'][2:] # skip scale, background
422                     if p[0] in set(fixed_1d + pd_1d)]
423    iqxy_parameters = [p[0]
424                       for p in info['parameters'][2:] # skip scale, background
425                       if p[0] in set(fixed_2d + pd_2d)]
426    volume_parameters = [p[0]
427                         for p in info['parameters']
428                         if p[4] == 'volume']
429
430    # Fill in defintions for volume parameters
431    if volume_parameters:
432        defines.append(('VOLUME_PARAMETERS',
433                        ','.join(volume_parameters)))
434        defines.append(('VOLUME_WEIGHT_PRODUCT',
435                        '*'.join(p + '_w' for p in volume_parameters)))
436
437    # Generate form_volume function from body only
438    if info['form_volume'] is not None:
439        if volume_parameters:
440            vol_par_decl = ', '.join('double ' + p for p in volume_parameters)
441        else:
442            vol_par_decl = 'void'
443        defines.append(('VOLUME_PARAMETER_DECLARATIONS',
444                        vol_par_decl))
445        fn = """\
446double form_volume(VOLUME_PARAMETER_DECLARATIONS);
447double form_volume(VOLUME_PARAMETER_DECLARATIONS) {
448    %(body)s
449}
450""" % {'body':info['form_volume']}
451        source.append(fn)
452
453    # Fill in definitions for Iq parameters
454    defines.append(('IQ_KERNEL_NAME', info['name'] + '_Iq'))
455    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters)))
456    if fixed_1d:
457        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS',
458                        ', \\\n    '.join('const double %s' % p for p in fixed_1d)))
459    if pd_1d:
460        defines.append(('IQ_WEIGHT_PRODUCT',
461                        '*'.join(p + '_w' for p in pd_1d)))
462        defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS',
463                        ', \\\n    '.join('const int N%s' % p for p in pd_1d)))
464        defines.append(('IQ_DISPERSION_LENGTH_SUM',
465                        '+'.join('N' + p for p in pd_1d)))
466        open_loops, close_loops = build_polydispersity_loops(pd_1d)
467        defines.append(('IQ_OPEN_LOOPS',
468                        open_loops.replace('\n', ' \\\n')))
469        defines.append(('IQ_CLOSE_LOOPS',
470                        close_loops.replace('\n', ' \\\n')))
471    if info['Iq'] is not None:
472        defines.append(('IQ_PARAMETER_DECLARATIONS',
473                        ', '.join('double ' + p for p in iq_parameters)))
474        fn = """\
475double Iq(double q, IQ_PARAMETER_DECLARATIONS);
476double Iq(double q, IQ_PARAMETER_DECLARATIONS) {
477    %(body)s
478}
479""" % {'body':info['Iq']}
480        source.append(fn)
481
482    # Fill in definitions for Iqxy parameters
483    defines.append(('IQXY_KERNEL_NAME', info['name'] + '_Iqxy'))
484    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters)))
485    if fixed_2d:
486        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS',
487                        ', \\\n    '.join('const double %s' % p for p in fixed_2d)))
488    if pd_2d:
489        defines.append(('IQXY_WEIGHT_PRODUCT',
490                        '*'.join(p + '_w' for p in pd_2d)))
491        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS',
492                        ', \\\n    '.join('const int N%s' % p for p in pd_2d)))
493        defines.append(('IQXY_DISPERSION_LENGTH_SUM',
494                        '+'.join('N' + p for p in pd_2d)))
495        open_loops, close_loops = build_polydispersity_loops(pd_2d)
496        defines.append(('IQXY_OPEN_LOOPS',
497                        open_loops.replace('\n', ' \\\n')))
498        defines.append(('IQXY_CLOSE_LOOPS',
499                        close_loops.replace('\n', ' \\\n')))
500    if info['Iqxy'] is not None:
501        defines.append(('IQXY_PARAMETER_DECLARATIONS',
502                        ', '.join('double ' + p for p in iqxy_parameters)))
503        fn = """\
504double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS);
505double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) {
506    %(body)s
507}
508""" % {'body':info['Iqxy']}
509        source.append(fn)
510
511    # Need to know if we have a theta parameter for Iqxy; it is not there
512    # for the magnetic sphere model, for example, which has a magnetic
513    # orientation but no shape orientation.
514    if 'theta' in pd_2d:
515        defines.append(('IQXY_HAS_THETA', '1'))
516
517    #for d in defines: print d
518    DEFINES = '\n'.join('#define %s %s' % (k, v) for k, v in defines)
519    SOURCES = '\n\n'.join(source)
520    return C_KERNEL_TEMPLATE % {
521        'DEFINES':DEFINES,
522        'SOURCES':SOURCES,
523        }
524
525def make(kernel_module):
526    """
527    Build an OpenCL/ctypes function from the definition in *kernel_module*.
528
529    The module can be loaded with a normal python import statement if you
530    know which module you need, or with __import__('sasmodels.model.'+name)
531    if the name is in a string.
532    """
533    # TODO: allow Iq and Iqxy to be defined in python
534    #print kernelfile
535    info = dict(
536        filename=abspath(kernel_module.__file__),
537        name=kernel_module.name,
538        title=kernel_module.title,
539        description=kernel_module.description,
540        parameters=COMMON_PARAMETERS + kernel_module.parameters,
541        source=getattr(kernel_module, 'source', []),
542        oldname=kernel_module.oldname,
543        oldpars=kernel_module.oldpars,
544        )
545    # Fill in attributes which default to None
546    info.update((k, getattr(kernel_module, k, None))
547                for k in ('ER', 'VR', 'form_volume', 'Iq', 'Iqxy'))
548    # Fill in the derived attributes
549    info['limits'] = dict((p[0], p[3]) for p in info['parameters'])
550    info['partype'] = categorize_parameters(info['parameters'])
551    info['defaults'] = dict((p[0], p[2]) for p in info['parameters'])
552
553    # Assume if one part of the kernel is python then all parts are.
554    source = make_model(info) if not callable(info['Iq']) else None
555    return source, info
556
557def doc(kernel_module):
558    """
559    Return the documentation for the model.
560    """
561    subst = dict(name=kernel_module.name.replace('_', '-'),
562                 label=" ".join(kernel_module.name.split('_')).capitalize(),
563                 title=kernel_module.title,
564                 parameters=make_partable(kernel_module.parameters),
565                 docs=kernel_module.__doc__)
566    return DOC_HEADER % subst
567
568
569
570def demo_time():
571    from .models import cylinder
572    import datetime
573    tic = datetime.datetime.now()
574    make(cylinder)
575    toc = (datetime.datetime.now() - tic).total_seconds()
576    print "time:", toc
577
578def main():
579    if len(sys.argv) <= 1:
580        print "usage: python -m sasmodels.generate modelname"
581    else:
582        name = sys.argv[1]
583        import sasmodels.models
584        __import__('sasmodels.models.' + name)
585        model = getattr(sasmodels.models, name)
586        source, _ = make(model)
587        print source
588
589if __name__ == "__main__":
590    main()
Note: See TracBrowser for help on using the repository browser.