source: sasmodels/sasmodels/generate.py @ 29fc2a3

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 29fc2a3 was 3c56da87, checked in by Paul Kienzle <pkienzle@…>, 9 years ago

lint cleanup

  • Property mode set to 100644
File size: 22.5 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 make_partable(pars):
250    """
251    Generate the parameter table to include in the sphinx documentation.
252    """
253    pars = COMMON_PARAMETERS + pars
254    column_widths = [
255        max(len(p[0]) for p in pars),
256        max(len(p[-1]) for p in pars),
257        max(len(RST_UNITS[p[1]]) for p in pars),
258        PARTABLE_VALUE_WIDTH,
259        ]
260    column_widths = [max(w, len(h))
261                     for w, h in zip(column_widths, PARTABLE_HEADERS)]
262
263    sep = " ".join("="*w for w in column_widths)
264    lines = [
265        sep,
266        " ".join("%-*s" % (w, h) for w, h in zip(column_widths, PARTABLE_HEADERS)),
267        sep,
268        ]
269    for p in pars:
270        lines.append(" ".join([
271            "%-*s" % (column_widths[0], p[0]),
272            "%-*s" % (column_widths[1], p[-1]),
273            "%-*s" % (column_widths[2], RST_UNITS[p[1]]),
274            "%*g" % (column_widths[3], p[2]),
275            ]))
276    lines.append(sep)
277    return "\n".join(lines)
278
279def _search(search_path, filename):
280    """
281    Find *filename* in *search_path*.
282
283    Raises ValueError if file does not exist.
284    """
285    for path in search_path:
286        target = joinpath(path, filename)
287        if exists(target):
288            return target
289    raise ValueError("%r not found in %s" % (filename, search_path))
290
291def sources(info):
292    """
293    Return a list of the sources file paths for the module.
294    """
295    search_path = [dirname(info['filename']),
296                   abspath(joinpath(dirname(__file__), 'models'))]
297    return [_search(search_path, f) for f in info['source']]
298
299def use_single(source):
300    """
301    Convert code from double precision to single precision.
302    """
303    # Convert double keyword to float.  Accept an 'n' parameter for vector
304    # values, where n is 2, 4, 8 or 16. Assume complex numbers are represented
305    # as cdouble which is typedef'd to double2.
306    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
307                    r'\1float\2', source)
308    # Convert floating point constants to single by adding 'f' to the end.
309    # OS/X driver complains if you don't do this.
310    source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?',
311                    r'\g<0>f', source)
312    return source
313
314
315def kernel_name(info, is_2D):
316    """
317    Name of the exported kernel symbol.
318    """
319    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq")
320
321
322def categorize_parameters(pars):
323    """
324    Build parameter categories out of the the parameter definitions.
325
326    Returns a dictionary of categories.
327    """
328    partype = {
329        'volume': [], 'orientation': [], 'magnetic': [], '': [],
330        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [],
331        'pd-rel': set(),
332    }
333
334    for p in pars:
335        name, ptype = p[0], p[4]
336        if ptype == 'volume':
337            partype['pd-1d'].append(name)
338            partype['pd-2d'].append(name)
339            partype['pd-rel'].add(name)
340        elif ptype == 'magnetic':
341            partype['fixed-2d'].append(name)
342        elif ptype == 'orientation':
343            partype['pd-2d'].append(name)
344        elif ptype == '':
345            partype['fixed-1d'].append(name)
346            partype['fixed-2d'].append(name)
347        else:
348            raise ValueError("unknown parameter type %r" % ptype)
349        partype[ptype].append(name)
350
351    return partype
352
353def indent(s, depth):
354    """
355    Indent a string of text with *depth* additional spaces on each line.
356    """
357    spaces = " "*depth
358    sep = "\n" + spaces
359    return spaces + sep.join(s.split("\n"))
360
361
362def build_polydispersity_loops(pd_pars):
363    """
364    Build polydispersity loops
365
366    Returns loop opening and loop closing
367    """
368    LOOP_OPEN = """\
369for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) {
370  const double %(name)s = loops[2*(%(name)s_i%(offset)s)];
371  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\
372"""
373    depth = 4
374    offset = ""
375    loop_head = []
376    loop_end = []
377    for name in pd_pars:
378        subst = {'name': name, 'offset': offset}
379        loop_head.append(indent(LOOP_OPEN % subst, depth))
380        loop_end.insert(0, (" "*depth) + "}")
381        offset += '+N' + name
382        depth += 2
383    return "\n".join(loop_head), "\n".join(loop_end)
384
385C_KERNEL_TEMPLATE = None
386def make_model(info):
387    """
388    Generate the code for the kernel defined by info, using source files
389    found in the given search path.
390    """
391    # TODO: need something other than volume to indicate dispersion parameters
392    # No volume normalization despite having a volume parameter.
393    # Thickness is labelled a volume in order to trigger polydispersity.
394    # May want a separate dispersion flag, or perhaps a separate category for
395    # disperse, but not volume.  Volume parameters also use relative values
396    # for the distribution rather than the absolute values used by angular
397    # dispersion.  Need to be careful that necessary parameters are available
398    # for computing volume even if we allow non-disperse volume parameters.
399
400    # Load template
401    global C_KERNEL_TEMPLATE
402    if C_KERNEL_TEMPLATE is None:
403        with open(C_KERNEL_TEMPLATE_PATH) as fid:
404            C_KERNEL_TEMPLATE = fid.read()
405
406    # Load additional sources
407    source = [open(f).read() for f in sources(info)]
408
409    # Prepare defines
410    defines = []
411    partype = info['partype']
412    pd_1d = partype['pd-1d']
413    pd_2d = partype['pd-2d']
414    fixed_1d = partype['fixed-1d']
415    fixed_2d = partype['fixed-1d']
416
417    iq_parameters = [p[0]
418                     for p in info['parameters'][2:] # skip scale, background
419                     if p[0] in set(fixed_1d + pd_1d)]
420    iqxy_parameters = [p[0]
421                       for p in info['parameters'][2:] # skip scale, background
422                       if p[0] in set(fixed_2d + pd_2d)]
423    volume_parameters = [p[0]
424                         for p in info['parameters']
425                         if p[4] == 'volume']
426
427    # Fill in defintions for volume parameters
428    if volume_parameters:
429        defines.append(('VOLUME_PARAMETERS',
430                        ','.join(volume_parameters)))
431        defines.append(('VOLUME_WEIGHT_PRODUCT',
432                        '*'.join(p + '_w' for p in volume_parameters)))
433
434    # Generate form_volume function from body only
435    if info['form_volume'] is not None:
436        if volume_parameters:
437            vol_par_decl = ', '.join('double ' + p for p in volume_parameters)
438        else:
439            vol_par_decl = 'void'
440        defines.append(('VOLUME_PARAMETER_DECLARATIONS',
441                        vol_par_decl))
442        fn = """\
443double form_volume(VOLUME_PARAMETER_DECLARATIONS);
444double form_volume(VOLUME_PARAMETER_DECLARATIONS) {
445    %(body)s
446}
447""" % {'body':info['form_volume']}
448        source.append(fn)
449
450    # Fill in definitions for Iq parameters
451    defines.append(('IQ_KERNEL_NAME', info['name'] + '_Iq'))
452    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters)))
453    if fixed_1d:
454        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS',
455                        ', \\\n    '.join('const double %s' % p for p in fixed_1d)))
456    if pd_1d:
457        defines.append(('IQ_WEIGHT_PRODUCT',
458                        '*'.join(p + '_w' for p in pd_1d)))
459        defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS',
460                        ', \\\n    '.join('const int N%s' % p for p in pd_1d)))
461        defines.append(('IQ_DISPERSION_LENGTH_SUM',
462                        '+'.join('N' + p for p in pd_1d)))
463        open_loops, close_loops = build_polydispersity_loops(pd_1d)
464        defines.append(('IQ_OPEN_LOOPS',
465                        open_loops.replace('\n', ' \\\n')))
466        defines.append(('IQ_CLOSE_LOOPS',
467                        close_loops.replace('\n', ' \\\n')))
468    if info['Iq'] is not None:
469        defines.append(('IQ_PARAMETER_DECLARATIONS',
470                        ', '.join('double ' + p for p in iq_parameters)))
471        fn = """\
472double Iq(double q, IQ_PARAMETER_DECLARATIONS);
473double Iq(double q, IQ_PARAMETER_DECLARATIONS) {
474    %(body)s
475}
476""" % {'body':info['Iq']}
477        source.append(fn)
478
479    # Fill in definitions for Iqxy parameters
480    defines.append(('IQXY_KERNEL_NAME', info['name'] + '_Iqxy'))
481    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters)))
482    if fixed_2d:
483        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS',
484                        ', \\\n    '.join('const double %s' % p for p in fixed_2d)))
485    if pd_2d:
486        defines.append(('IQXY_WEIGHT_PRODUCT',
487                        '*'.join(p + '_w' for p in pd_2d)))
488        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS',
489                        ', \\\n    '.join('const int N%s' % p for p in pd_2d)))
490        defines.append(('IQXY_DISPERSION_LENGTH_SUM',
491                        '+'.join('N' + p for p in pd_2d)))
492        open_loops, close_loops = build_polydispersity_loops(pd_2d)
493        defines.append(('IQXY_OPEN_LOOPS',
494                        open_loops.replace('\n', ' \\\n')))
495        defines.append(('IQXY_CLOSE_LOOPS',
496                        close_loops.replace('\n', ' \\\n')))
497    if info['Iqxy'] is not None:
498        defines.append(('IQXY_PARAMETER_DECLARATIONS',
499                        ', '.join('double ' + p for p in iqxy_parameters)))
500        fn = """\
501double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS);
502double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) {
503    %(body)s
504}
505""" % {'body':info['Iqxy']}
506        source.append(fn)
507
508    # Need to know if we have a theta parameter for Iqxy; it is not there
509    # for the magnetic sphere model, for example, which has a magnetic
510    # orientation but no shape orientation.
511    if 'theta' in pd_2d:
512        defines.append(('IQXY_HAS_THETA', '1'))
513
514    #for d in defines: print d
515    DEFINES = '\n'.join('#define %s %s' % (k, v) for k, v in defines)
516    SOURCES = '\n\n'.join(source)
517    return C_KERNEL_TEMPLATE % {
518        'DEFINES':DEFINES,
519        'SOURCES':SOURCES,
520        }
521
522def make(kernel_module):
523    """
524    Build an OpenCL/ctypes function from the definition in *kernel_module*.
525
526    The module can be loaded with a normal python import statement if you
527    know which module you need, or with __import__('sasmodels.model.'+name)
528    if the name is in a string.
529    """
530    # TODO: allow Iq and Iqxy to be defined in python
531    #print kernelfile
532    info = dict(
533        filename=abspath(kernel_module.__file__),
534        name=kernel_module.name,
535        title=kernel_module.title,
536        description=kernel_module.description,
537        parameters=COMMON_PARAMETERS + kernel_module.parameters,
538        source=getattr(kernel_module, 'source', []),
539        oldname=kernel_module.oldname,
540        oldpars=kernel_module.oldpars,
541        )
542    # Fill in attributes which default to None
543    info.update((k, getattr(kernel_module, k, None))
544                for k in ('ER', 'VR', 'form_volume', 'Iq', 'Iqxy'))
545    # Fill in the derived attributes
546    info['limits'] = dict((p[0], p[3]) for p in info['parameters'])
547    info['partype'] = categorize_parameters(info['parameters'])
548    info['defaults'] = dict((p[0], p[2]) for p in info['parameters'])
549
550    # Assume if one part of the kernel is python then all parts are.
551    source = make_model(info) if not callable(info['Iq']) else None
552    return source, info
553
554def doc(kernel_module):
555    """
556    Return the documentation for the model.
557    """
558    subst = dict(name=kernel_module.name.replace('_', '-'),
559                 label=" ".join(kernel_module.name.split('_')).capitalize(),
560                 title=kernel_module.title,
561                 parameters=make_partable(kernel_module.parameters),
562                 docs=kernel_module.__doc__)
563    return DOC_HEADER % subst
564
565
566
567def demo_time():
568    from .models import cylinder
569    import datetime
570    tic = datetime.datetime.now()
571    make(cylinder)
572    toc = (datetime.datetime.now() - tic).total_seconds()
573    print "time:", toc
574
575def main():
576    if len(sys.argv) <= 1:
577        print "usage: python -m sasmodels.generate modelname"
578    else:
579        name = sys.argv[1]
580        import sasmodels.models
581        __import__('sasmodels.models.' + name)
582        model = getattr(sasmodels.models, name)
583        source, _ = make(model)
584        print source
585
586if __name__ == "__main__":
587    main()
Note: See TracBrowser for help on using the repository browser.