source: sasmodels/sasmodels/generate.py @ a8a1d48

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since a8a1d48 was a8a1d48, checked in by GitHub <noreply@…>, 4 months ago

Set structure factor background default to 0.0. Closes #1157.

  • Property mode set to 100644
File size: 39.7 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    *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc
10    for a rotationally symmetric form with particular dimensions.
11    qab, qc are determined from shape orientation and scattering angles.
12    This call is used if the shape has orientation parameters theta and phi.
13
14    *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc
15    for a form with particular dimensions.  qa, qb, qc are determined from
16    shape orientation and scattering angles. This call is used if the shape
17    has orientation parameters theta, phi and psi.
18
19    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy.  Use this
20    to create an arbitrary 2D theory function, needed for q-dependent
21    background functions and for models with non-uniform magnetism.
22
23    *form_volume(p1, p2, ...)* returns the volume of the form with particular
24    dimension, or 1.0 if no volume normalization is required.
25
26    *shell_volume(p1, p2, ...)* returns the volume of the shell for forms
27    which are hollow.
28
29    *effective_radius(mode, p1, p2, ...)* returns the effective radius of
30    the form with particular dimensions.  Mode determines the type of
31    effective radius returned, with mode=1 for equivalent volume.
32
33    #define INVALID(v) (expr)  returns False if v.parameter is invalid
34    for some parameter or other (e.g., v.bell_radius < v.radius).  If
35    necessary, the expression can call a function.
36
37These functions are defined in a kernel module .py script and an associated
38set of .c files.  The model constructor will use them to create models with
39polydispersity across volume and orientation parameters, and provide
40scale and background parameters for each model.
41
42C code should be stylized C-99 functions written for OpenCL. All functions
43need prototype declarations even if the are defined before they are used.
44Although OpenCL supports *#include* preprocessor directives, the list of
45includes should be given as part of the metadata in the kernel module
46definition. The included files should be listed using a path relative to the
47kernel module, or if using "lib/file.c" if it is one of the standard includes
48provided with the sasmodels source. The includes need to be listed in order
49so that functions are defined before they are used.
50
51Floating point values should be declared as *double*.  For single precision
52calculations, *double* will be replaced by *float*.  The single precision
53conversion will also tag floating point constants with "f" to make them
54single precision constants.  When using integral values in floating point
55expressions, they should be expressed as floating point values by including
56a decimal point.  This includes 0., 1. and 2.
57
58OpenCL has a *sincos* function which can improve performance when both
59the *sin* and *cos* values are needed for a particular argument.  Since
60this function does not exist in C99, all use of *sincos* should be
61replaced by the macro *SINCOS(value, sn, cn)* where *sn* and *cn* are
62previously declared *double* variables.  When compiled for systems without
63OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.   If *value* is
64an expression, it will appear twice in this case; whether or not it will be
65evaluated twice depends on the quality of the compiler.
66
67If the input parameters are invalid, the scattering calculator should
68return a negative number. Particularly with polydispersity, there are
69some sets of shape parameters which lead to nonsensical forms, such
70as a capped cylinder where the cap radius is smaller than the
71cylinder radius.  The polydispersity calculation will ignore these points,
72effectively chopping the parameter weight distributions at the boundary
73of the infeasible region.  The resulting scattering will be set to
74background.  This will work correctly even when polydispersity is off.
75
76The kernel module must set variables defining the kernel meta data:
77
78    *id* is an implicit variable formed from the filename.  It will be
79    a valid python identifier, and will be used as the reference into
80    the html documentation, with '_' replaced by '-'.
81
82    *name* is the model name as displayed to the user.  If it is missing,
83    it will be constructed from the id.
84
85    *title* is a short description of the model, suitable for a tool tip,
86    or a one line model summary in a table of models.
87
88    *description* is an extended description of the model to be displayed
89    while the model parameters are being edited.
90
91    *parameters* is the list of parameters.  Parameters in the kernel
92    functions must appear in the same order as they appear in the
93    parameters list.  Two additional parameters, *scale* and *background*
94    are added to the beginning of the parameter list.  They will show up
95    in the documentation as model parameters, but they are never sent to
96    the kernel functions.  Note that *effect_radius* and *volfraction*
97    must occur first in structure factor calculations.
98
99    *category* is the default category for the model.  The category is
100    two level structure, with the form "group:section", indicating where
101    in the manual the model will be located.  Models are alphabetical
102    within their section.
103
104    *source* is the list of C-99 source files that must be joined to
105    create the OpenCL kernel functions.  The files defining the functions
106    need to be listed before the files which use the functions.
107
108    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing
109    the C source code for the body of the volume, Iq, and Iqac functions
110    respectively.  These can also be defined in the last source file.
111
112    *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the
113    kernel.  If they are marked as *Iq.vectorized = True* then the
114    kernel is passed the entire *q* vector at once, otherwise it is
115    passed values one *q* at a time.  The performance improvement of
116    this step is significant.
117
118    *demo* is a dictionary of parameter=value defining a set of
119    parameters to use by default when *compare* is called.  Any
120    parameter not set in *demo* gets the initial value from the
121    parameter list.  *demo* is mostly needed to set the default
122    polydispersity values for tests.
123
124A :class:`modelinfo.ModelInfo` structure is constructed from the kernel meta
125data and returned to the caller.
126
127The doc string at the start of the kernel module will be used to
128construct the model documentation web pages.  Embedded figures should
129appear in the subdirectory "img" beside the model definition, and tagged
130with the kernel module name to avoid collision with other models.  Some
131file systems are case-sensitive, so only use lower case characters for
132file names and extensions.
133
134Code follows the C99 standard with the following extensions and conditions::
135
136    M_PI_180 = pi/180
137    M_4PI_3 = 4pi/3
138    square(x) = x*x
139    cube(x) = x*x*x
140    sas_sinx_x(x) = sin(x)/x, with sin(0)/0 -> 1
141    all double precision constants must include the decimal point
142    all double declarations may be converted to half, float, or long double
143    FLOAT_SIZE is the number of bytes in the converted variables
144
145:func:`load_kernel_module` loads the model definition file and
146:func:`modelinfo.make_model_info` parses it. :func:`make_source`
147converts C-based model definitions to C source code, including the
148polydispersity integral.  :func:`model_sources` returns the list of
149source files the model depends on, and :func:`timestamp` returns
150the latest time stamp amongst the source files (so you can check if
151the model needs to be rebuilt).
152
153The function :func:`make_doc` extracts the doc string and adds the
154parameter table to the top.  *make_figure* in *sasmodels/doc/genmodel*
155creates the default figure for the model.  [These two sets of code
156should mignrate into docs.py so docs can be updated in one place].
157"""
158from __future__ import print_function
159
160# TODO: determine which functions are useful outside of generate
161#__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
162
163import sys
164from os import environ
165from os.path import abspath, dirname, join as joinpath, exists, getmtime, sep
166import re
167import string
168from zlib import crc32
169from inspect import currentframe, getframeinfo
170import logging
171
172import numpy as np  # type: ignore
173
174from .modelinfo import Parameter
175from .custom import load_custom_kernel_module
176
177# pylint: disable=unused-import
178try:
179    from typing import Tuple, Sequence, Iterator, Dict
180    from .modelinfo import ModelInfo
181except ImportError:
182    pass
183# pylint: enable=unused-import
184
185logger = logging.getLogger(__name__)
186
187# jitter projection to use in the kernel code.  See explore/jitter.py
188# for details.  To change it from a program, set generate.PROJECTION.
189PROJECTION = 1
190
191def get_data_path(external_dir, target_file):
192    path = abspath(dirname(__file__))
193    if exists(joinpath(path, target_file)):
194        return path
195
196    # check next to exe/zip file
197    exepath = dirname(sys.executable)
198    path = joinpath(exepath, external_dir)
199    if exists(joinpath(path, target_file)):
200        return path
201
202    # check in py2app Contents/Resources
203    path = joinpath(exepath, '..', 'Resources', external_dir)
204    if exists(joinpath(path, target_file)):
205        return abspath(path)
206
207    raise RuntimeError('Could not find '+joinpath(external_dir, target_file))
208
209EXTERNAL_DIR = 'sasmodels-data'
210DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_iq.c')
211MODEL_PATH = joinpath(DATA_PATH, 'models')
212
213F16 = np.dtype('float16')
214F32 = np.dtype('float32')
215F64 = np.dtype('float64')
216try:  # CRUFT: older numpy does not support float128
217    F128 = np.dtype('float128')
218except TypeError:
219    F128 = None
220
221# Conversion from units defined in the parameter table for each model
222# to units displayed in the sphinx documentation.
223# This section associates the unit with the macro to use to produce the LaTex
224# code.  The macro itself needs to be defined in sasmodels/doc/rst_prolog.
225#
226# NOTE: there is an RST_PROLOG at the end of this file which is NOT
227# used for the bundled documentation. Still as long as we are defining the macros
228# in two places any new addition should define the macro in both places.
229RST_UNITS = {
230    "Ang": "|Ang|",
231    "1/Ang": "|Ang^-1|",
232    "1/Ang^2": "|Ang^-2|",
233    "Ang^3": "|Ang^3|",
234    "Ang^2": "|Ang^2|",
235    "1e15/cm^3": "|1e15cm^3|",
236    "Ang^3/mol": "|Ang^3|/mol",
237    "1e-6/Ang^2": "|1e-6Ang^-2|",
238    "degrees": "degree",
239    "1/cm": "|cm^-1|",
240    "Ang/cm": "|Ang*cm^-1|",
241    "g/cm^3": "|g/cm^3|",
242    "mg/m^2": "|mg/m^2|",
243    "": "None",
244    }
245
246# Headers for the parameters tables in th sphinx documentation
247PARTABLE_HEADERS = [
248    "Parameter",
249    "Description",
250    "Units",
251    "Default value",
252    ]
253
254# Minimum width for a default value (this is shorter than the column header
255# width, so will be ignored).
256PARTABLE_VALUE_WIDTH = 10
257
258# Documentation header for the module, giving the model name, its short
259# description and its parameter table.  The remainder of the doc comes
260# from the module docstring.
261DOC_HEADER = """.. _%(id)s:
262
263%(name)s
264=======================================================
265
266%(title)s
267
268%(parameters)s
269
270%(returns)s
271
272%(docs)s
273"""
274
275
276def set_integration_size(info, n):
277    # type: (ModelInfo, int) -> None
278    """
279    Update the model definition, replacing the gaussian integration with
280    a gaussian integration of a different size.
281
282    Note: this really ought to be a method in modelinfo, but that leads to
283    import loops.
284    """
285    if info.source and any(lib.startswith('lib/gauss') for lib in info.source):
286        from .gengauss import gengauss
287        path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n)
288        if not exists(path):
289            gengauss(n, path)
290        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')
291                       else lib for lib in info.source]
292
293def format_units(units):
294    # type: (str) -> str
295    """
296    Convert units into ReStructured Text format.
297    """
298    return "string" if isinstance(units, list) else RST_UNITS.get(units, units)
299
300
301def make_partable(pars):
302    # type: (List[Parameter]) -> str
303    """
304    Generate the parameter table to include in the sphinx documentation.
305    """
306    column_widths = [
307        max(len(p.name) for p in pars),
308        max(len(p.description) for p in pars),
309        max(len(format_units(p.units)) for p in pars),
310        PARTABLE_VALUE_WIDTH,
311        ]
312    column_widths = [max(w, len(h))
313                     for w, h in zip(column_widths, PARTABLE_HEADERS)]
314
315    underbar = " ".join("="*w for w in column_widths)
316    lines = [
317        underbar,
318        " ".join("%-*s" % (w, h)
319                 for w, h in zip(column_widths, PARTABLE_HEADERS)),
320        underbar,
321        ]
322    for p in pars:
323        lines.append(" ".join([
324            "%-*s" % (column_widths[0], p.name),
325            "%-*s" % (column_widths[1], p.description),
326            "%-*s" % (column_widths[2], format_units(p.units)),
327            "%*g" % (column_widths[3], p.default),
328            ]))
329    lines.append(underbar)
330    return "\n".join(lines)
331
332
333def _search(search_path, filename):
334    # type: (List[str], str) -> str
335    """
336    Find *filename* in *search_path*.
337
338    Raises ValueError if file does not exist.
339    """
340    for path in search_path:
341        target = joinpath(path, filename)
342        if exists(target):
343            return target
344    raise ValueError("%r not found in %s" % (filename, search_path))
345
346
347def model_sources(model_info):
348    # type: (ModelInfo) -> List[str]
349    """
350    Return a list of the sources file paths for the module.
351    """
352    search_path = [dirname(model_info.filename), MODEL_PATH]
353    return [_search(search_path, f) for f in model_info.source]
354
355
356def dll_timestamp(model_info):
357    # type: (ModelInfo) -> int
358    """
359    Return a timestamp for the model corresponding to the most recently
360    changed file or dependency.
361    """
362    # TODO: fails DRY; templates appear two places.
363    model_templates = [joinpath(DATA_PATH, filename)
364                       for filename in ('kernel_header.c', 'kernel_iq.c')]
365    source_files = (model_sources(model_info)
366                    + model_templates
367                    + [model_info.filename])
368    # Note: file may not exist when it is a standard model from library.zip
369    times = [getmtime(f) for f in source_files if exists(f)]
370    newest = max(times) if times else 0
371    return newest
372
373def ocl_timestamp(model_info):
374    # type: (ModelInfo) -> int
375    """
376    Return a timestamp for the model corresponding to the most recently
377    changed file or dependency.
378
379    Note that this does not look at the time stamps for the OpenCL header
380    information since that need not trigger a recompile of the DLL.
381    """
382    # TODO: fails DRY; templates appear two places.
383    model_templates = [joinpath(DATA_PATH, filename)
384                       for filename in ('kernel_header.c', 'kernel_iq.c')]
385    source_files = (model_sources(model_info)
386                    + model_templates
387                    + [model_info.filename])
388    # Note: file may not exist when it is a standard model from library.zip
389    times = [getmtime(f) for f in source_files if exists(f)]
390    newest = max(times) if times else 0
391    return newest
392
393def tag_source(source):
394    # type: (str) -> str
395    """
396    Return a unique tag for the source code.
397    """
398    # Note: need 0xffffffff&val to force an unsigned 32-bit number
399    try:
400        source = source.encode('utf8')
401    except AttributeError: # bytes has no encode attribute in python 3
402        pass
403    return "%08X"%(0xffffffff&crc32(source))
404
405def convert_type(source, dtype):
406    # type: (str, np.dtype) -> str
407    """
408    Convert code from double precision to the desired type.
409
410    Floating point constants are tagged with 'f' for single precision or 'L'
411    for long double precision.
412    """
413    source = _fix_tgmath_int(source)
414    if dtype == F16:
415        fbytes = 2
416        source = _convert_type(source, "half", "f")
417    elif dtype == F32:
418        fbytes = 4
419        source = _convert_type(source, "float", "f")
420    elif dtype == F64:
421        fbytes = 8
422        # no need to convert if it is already double
423    elif dtype == F128:
424        fbytes = 16
425        source = _convert_type(source, "long double", "L")
426    else:
427        raise ValueError("Unexpected dtype in source conversion: %s" % dtype)
428    return ("#define FLOAT_SIZE %d\n" % fbytes)+source
429
430
431def _convert_type(source, type_name, constant_flag):
432    # type: (str, str, str) -> str
433    """
434    Replace 'double' with *type_name* in *source*, tagging floating point
435    constants with *constant_flag*.
436    """
437    # Convert double keyword to float/long double/half.
438    # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16.
439    # Assume complex numbers are represented as cdouble which is typedef'd
440    # to double2.
441    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
442                    r'\1%s\2'%type_name, source)
443    source = _tag_float(source, constant_flag)
444    return source
445
446TGMATH_INT_RE = re.compile(r"""
447(?: # Non-capturing match; not lookbehind since pattern length is variable
448  \b              # word boundary
449   # various math functions
450  (a?(sin|cos|tan)h? | atan2
451   | erfc? | tgamma
452   | exp(2|10|m1)? | log(2|10|1p)? | pow[nr]? | sqrt | rsqrt | rootn
453   | fabs | fmax | fmin
454   )
455  \s*[(]\s*       # open parenthesis
456)
457[+-]?(0|[1-9]\d*) # integer
458(?=               # lookahead match: don't want to move from end of int
459  \s*[,)]         # comma or close parenthesis for end of argument
460)                 # end lookahead
461""", re.VERBOSE)
462def _fix_tgmath_int(source):
463    # type: (str) -> str
464    """
465    Replace f(integer) with f(integer.) for sin, cos, pow, etc.
466
467    OS X OpenCL complains that it can't resolve the type generic calls to
468    the standard math functions when they are called with integer constants,
469    but this does not happen with the Windows Intel driver for example.
470    To avoid confusion on the matrix marketplace, automatically promote
471    integers to floats if we recognize them in the source.
472
473    The specific functions we look for are:
474
475        trigonometric: sin, asin, sinh, asinh, etc., and atan2
476        exponential:   exp, exp2, exp10, expm1, log, log2, log10, logp1
477        power:         pow, pown, powr, sqrt, rsqrt, rootn
478        special:       erf, erfc, tgamma
479        float:         fabs, fmin, fmax
480
481    Note that we don't convert the second argument of dual argument
482    functions: atan2, fmax, fmin, pow, powr.  This could potentially
483    be a problem for pow(x, 2), but that case seems to work without change.
484    """
485    out = TGMATH_INT_RE.sub(r'\g<0>.', source)
486    return out
487
488
489# Floating point regular expression
490#
491# Define parts:
492#
493#    E = [eE][+-]?\d+    : Exponent
494#    P = [.]             : Decimal separator
495#    N = [1-9]\d*        : Natural number, no leading zeros
496#    Z = 0               : Zero
497#    F = \d+             : Fractional number, maybe leading zeros
498#    F? = \d*            : Optional fractional number
499#
500# We want to reject bare natural numbers and bare decimal points, so we
501# need to tediously outline the cases where we have either a fraction or
502# an exponent:
503#
504#   ( ZP | ZPF | ZE | ZPE | ZPFE | NP | NPF | NE | NPE | NPFE | PF | PFE )
505#
506#
507# We can then join cases by making parts optional.  The following are
508# some ways to do this:
509#
510#   ( (Z|N)(P|PF|E|PE|PFE) | PFE? )                   # Split on lead
511#     => ( (Z|N)(PF?|(PF?)?E) | PFE? )
512#   ( ((Z|N)PF?|PF)E? | (Z|N)E)                       # Split on point
513#   ( (ZP|ZPF|NP|NPF|PF) | (Z|ZP|ZPF|N|NP|NPF|PF)E )  # Split on E
514#     => ( ((Z|N)PF?|PF) | ((Z|N)(PF?)? | PF) E )
515FLOAT_RE = re.compile(r"""
516    (?<!\w)  # use negative lookbehind since '.' confuses \b test
517    # use split on lead to match float ( (Z|N)(PF?|(PF?)?E) | PFE? )
518    ( ( 0 | [1-9]\d* )                     # ( ( Z | N )
519      ([.]\d* | ([.]\d*)? [eE][+-]?\d+ )   #   (PF? | (PF?)? E )
520    | [.]\d+ ([eE][+-]?\d+)?               # | PF (E)?
521    )                                      # )
522    (?!\w)  # use negative lookahead since '.' confuses \b test
523    """, re.VERBOSE)
524def _tag_float(source, constant_flag):
525    # Convert floating point constants to single by adding 'f' to the end,
526    # or long double with an 'L' suffix.  OS/X complains if you don't do this.
527    out = FLOAT_RE.sub(r'\g<0>%s'%constant_flag, source)
528    #print("in",repr(source),"out",repr(out), constant_flag)
529    return out
530
531def test_tag_float():
532    """check that floating point constants are properly identified and tagged with 'f'"""
533
534    cases = """
535ZP  : 0.
536ZPF : 0.0,0.01,0.1
537Z  E: 0e+001
538ZP E: 0.E0
539ZPFE: 0.13e-031
540NP  : 1., 12.
541NPF : 1.0001, 1.1, 1.0
542N  E: 1e0, 37E-080
543NP E: 1.e0, 37.E-080
544NPFE: 845.017e+22
545 PF : .1, .0, .0100
546 PFE: .6e+9, .82E-004
547# isolated cases
5480.
5491e0
5500.13e-013
551# untouched
552struct3.e3, 03.05.67, 37
553# expressions
5543.75+-1.6e-7-27+13.2
555a3.e2 - 0.
5564*atan(1)
5574.*atan(1.)
558"""
559
560    output = """
561ZP  : 0.f
562ZPF : 0.0f,0.01f,0.1f
563Z  E: 0e+001f
564ZP E: 0.E0f
565ZPFE: 0.13e-031f
566NP  : 1.f, 12.f
567NPF : 1.0001f, 1.1f, 1.0f
568N  E: 1e0f, 37E-080f
569NP E: 1.e0f, 37.E-080f
570NPFE: 845.017e+22f
571 PF : .1f, .0f, .0100f
572 PFE: .6e+9f, .82E-004f
573# isolated cases
5740.f
5751e0f
5760.13e-013f
577# untouched
578struct3.e3, 03.05.67, 37
579# expressions
5803.75f+-1.6e-7f-27+13.2f
581a3.e2 - 0.f
5824*atan(1)
5834.f*atan(1.f)
584"""
585
586    for case_in, case_out in zip(cases.split('\n'), output.split('\n')):
587        out = _tag_float(case_in, 'f')
588        assert case_out == out, "%r => %r"%(case_in, out)
589
590
591def kernel_name(model_info, variant):
592    # type: (ModelInfo, str) -> str
593    """
594    Name of the exported kernel symbol.
595
596    *variant* is "Iq", "Iqxy" or "Imagnetic".
597    """
598    return model_info.name + "_" + variant
599
600
601def indent(s, depth):
602    # type: (str, int) -> str
603    """
604    Indent a string of text with *depth* additional spaces on each line.
605    """
606    spaces = " "*depth
607    interline_separator = "\n" + spaces
608    return spaces + interline_separator.join(s.split("\n"))
609
610
611_template_cache = {}  # type: Dict[str, Tuple[int, str, str]]
612def load_template(filename):
613    # type: (str) -> str
614    """
615    Load template file from sasmodels resource directory.
616    """
617    path = joinpath(DATA_PATH, filename)
618    mtime = getmtime(path)
619    if filename not in _template_cache or mtime > _template_cache[filename][0]:
620        with open(path) as fid:
621            _template_cache[filename] = (mtime, fid.read(), path)
622    return _template_cache[filename][1], path
623
624
625_FN_TEMPLATE = """\
626double %(name)s(%(pars)s);
627double %(name)s(%(pars)s) {
628#line %(line)d "%(filename)s"
629    %(body)s
630}
631
632"""
633def _gen_fn(model_info, name, pars):
634    # type: (ModelInfo, str, List[Parameter]) -> str
635    """
636    Generate a function given pars and body.
637
638    Returns the following string::
639
640         double fn(double a, double b, ...);
641         double fn(double a, double b, ...) {
642             ....
643         }
644    """
645    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
646    body = getattr(model_info, name)
647    filename = model_info.filename
648    # Note: if symbol is defined strangely in the module then default it to 1
649    lineno = model_info.lineno.get(name, 1)
650    return _FN_TEMPLATE % {
651        'name': name, 'pars': par_decl, 'body': body,
652        'filename': filename.replace('\\', '\\\\'), 'line': lineno,
653    }
654
655
656def _call_pars(prefix, pars):
657    # type: (str, List[Parameter]) -> List[str]
658    """
659    Return a list of *prefix+parameter* from parameter items.
660
661    *prefix* should be "v." if v is a struct.
662    """
663    return [p.as_call_reference(prefix) for p in pars]
664
665
666_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]",
667                           flags=re.MULTILINE)
668def find_xy_mode(source):
669    # type: (List[str]) -> bool
670    """
671    Return the xy mode as qa, qac, qabc or qxy.
672
673    Note this is not a C parser, and so can be easily confused by
674    non-standard syntax.  Also, it will incorrectly identify the following
675    as having 2D models::
676
677        /*
678        double Iqac(qab, qc, ...) { ... fill this in later ... }
679        */
680
681    If you want to comment out the function, use // on the front of the
682    line::
683
684        /*
685        // double Iqac(qab, qc, ...) { ... fill this in later ... }
686        */
687
688    """
689    for code in source:
690        m = _IQXY_PATTERN.search(code)
691        if m is not None:
692            return m.group('mode')
693    return 'qa'
694
695# Note: not presently used.  Need to know whether Fq is available before
696# trying to compile the source.  Leave the code here in case we decide that
697# define have_Fq for each form factor is too tedious and error prone.
698_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE)
699def contains_Fq(source):
700    # type: (List[str]) -> bool
701    """
702    Return True if C source defines "void Fq(".
703    """
704    for code in source:
705        if _FQ_PATTERN.search(code) is not None:
706            return True
707    return False
708
709_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE)
710def contains_shell_volume(source):
711    # type: (List[str]) -> bool
712    """
713    Return True if C source defines "double shell_volume(".
714    """
715    for code in source:
716        if _SHELL_VOLUME_PATTERN.search(code) is not None:
717            return True
718    return False
719
720def _add_source(source, code, path, lineno=1):
721    """
722    Add a file to the list of source code chunks, tagged with path and line.
723    """
724    path = path.replace('\\', '\\\\')
725    source.append('#line %d "%s"' % (lineno, path))
726    source.append(code)
727
728def make_source(model_info):
729    # type: (ModelInfo) -> Dict[str, str]
730    """
731    Generate the OpenCL/ctypes kernel from the module info.
732
733    Uses source files found in the given search path.  Returns None if this
734    is a pure python model, with no C source components.
735    """
736    if callable(model_info.Iq):
737        raise ValueError("can't compile python model")
738        #return None
739
740    # TODO: need something other than volume to indicate dispersion parameters
741    # No volume normalization despite having a volume parameter.
742    # Thickness is labelled a volume in order to trigger polydispersity.
743    # May want a separate dispersion flag, or perhaps a separate category for
744    # disperse, but not volume.  Volume parameters also use relative values
745    # for the distribution rather than the absolute values used by angular
746    # dispersion.  Need to be careful that necessary parameters are available
747    # for computing volume even if we allow non-disperse volume parameters.
748    partable = model_info.parameters
749
750    # Load templates and user code
751    kernel_header = load_template('kernel_header.c')
752    kernel_code = load_template('kernel_iq.c')
753    user_code = [(f, open(f).read()) for f in model_sources(model_info)]
754
755    # Build initial sources
756    source = []
757    _add_source(source, *kernel_header)
758    for path, code in user_code:
759        _add_source(source, code, path)
760    if model_info.c_code:
761        _add_source(source, model_info.c_code, model_info.filename,
762                    lineno=model_info.lineno.get('c_code', 1))
763
764    # Make parameters for q, qx, qy so that we can use them in declarations
765    q, qx, qy, qab, qa, qb, qc \
766        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
767
768    # Generate form_volume function, etc. from body only
769    if isinstance(model_info.form_volume, str):
770        pars = partable.form_volume_parameters
771        source.append(_gen_fn(model_info, 'form_volume', pars))
772    if isinstance(model_info.shell_volume, str):
773        pars = partable.form_volume_parameters
774        source.append(_gen_fn(model_info, 'shell_volume', pars))
775    if isinstance(model_info.Iq, str):
776        pars = [q] + partable.iq_parameters
777        source.append(_gen_fn(model_info, 'Iq', pars))
778    if isinstance(model_info.Iqxy, str):
779        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters
780        source.append(_gen_fn(model_info, 'Iqxy', pars))
781    if isinstance(model_info.Iqac, str):
782        pars = [qab, qc] + partable.iq_parameters
783        source.append(_gen_fn(model_info, 'Iqac', pars))
784    if isinstance(model_info.Iqabc, str):
785        pars = [qa, qb, qc] + partable.iq_parameters
786        source.append(_gen_fn(model_info, 'Iqabc', pars))
787
788    # Check for shell_volume in source
789    is_hollow = contains_shell_volume(source)
790
791    # What kind of 2D model do we need?  Is it consistent with the parameters?
792    xy_mode = find_xy_mode(source)
793    if xy_mode == 'qabc' and not partable.is_asymmetric:
794        raise ValueError("asymmetric oriented models need to define Iqabc")
795    elif xy_mode == 'qac' and partable.is_asymmetric:
796        raise ValueError("symmetric oriented models need to define Iqac")
797    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'):
798        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode)
799    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'):
800        if xy_mode == 'qxy':
801            logger.warn("oriented shapes should define Iqac or Iqabc")
802        else:
803            raise ValueError("Expected function Iqac or Iqabc for oriented shape")
804
805    # Define the parameter table
806    lineno = getframeinfo(currentframe()).lineno + 2
807    source.append('#line %d "sasmodels/generate.py"'%lineno)
808    #source.append('introduce breakage in generate to test lineno reporting')
809    source.append("#define PARAMETER_TABLE \\")
810    source.append("\\\n".join(p.as_definition()
811                              for p in partable.kernel_parameters))
812    # Define the function calls
813    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0"
814    if partable.form_volume_parameters:
815        refs = _call_pars("_v.", partable.form_volume_parameters)
816        if is_hollow:
817            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2)
818        else:
819            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs))
820        if model_info.effective_radius_type:
821            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs))
822    else:
823        # Model doesn't have volume.  We could make the kernel run a little
824        # faster by not using/transferring the volume normalizations, but
825        # the ifdef's reduce readability more than is worthwhile.
826        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)"
827    source.append(call_volume)
828    source.append(call_effective_radius)
829    model_refs = _call_pars("_v.", partable.iq_parameters)
830
831    if model_info.have_Fq:
832        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
833        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars
834        clear_iq = "#undef CALL_FQ"
835    else:
836        pars = ",".join(["_q"] + model_refs)
837        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
838        clear_iq = "#undef CALL_IQ"
839    if xy_mode == 'qabc':
840        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
841        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars
842        clear_iqxy = "#undef CALL_IQ_ABC"
843    elif xy_mode == 'qac':
844        pars = ",".join(["_qa", "_qc"] + model_refs)
845        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
846        clear_iqxy = "#undef CALL_IQ_AC"
847    elif xy_mode == 'qa' and not model_info.have_Fq:
848        pars = ",".join(["_qa"] + model_refs)
849        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
850        clear_iqxy = "#undef CALL_IQ_A"
851    elif xy_mode == 'qa' and model_info.have_Fq:
852        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs)
853        # Note: uses rare C construction (expr1, expr2) which computes
854        # expr1 then expr2 and evaluates to expr2.  This allows us to
855        # leave it looking like a function even though it is returning
856        # its values by reference.
857        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars
858        clear_iqxy = "#undef CALL_FQ_A"
859    elif xy_mode == 'qxy':
860        orientation_refs = _call_pars("_v.", partable.orientation_parameters)
861        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs)
862        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars
863        clear_iqxy = "#undef CALL_IQ_XY"
864        if partable.orientation_parameters:
865            call_iqxy += "\n#define HAVE_THETA"
866            clear_iqxy += "\n#undef HAVE_THETA"
867        if partable.is_asymmetric:
868            call_iqxy += "\n#define HAVE_PSI"
869            clear_iqxy += "\n#undef HAVE_PSI"
870
871
872    magpars = [k-2 for k, p in enumerate(partable.call_parameters)
873               if p.type == 'sld']
874    # Fill in definitions for numbers of parameters
875    source.append("#define MAX_PD %s"%partable.max_pd)
876    source.append("#define NUM_PARS %d"%partable.npars)
877    source.append("#define NUM_VALUES %d" % partable.nvalues)
878    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
879    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
880    source.append("#define PROJECTION %d"%PROJECTION)
881    # TODO: allow mixed python/opencl kernels?
882    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
883    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
884
885    result = {
886        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
887        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
888    }
889    return result
890
891
892def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name):
893    # type: ([str,str], str, str, str) -> List[str]
894    code = kernel[0]
895    path = kernel[1].replace('\\', '\\\\')
896    iq = [
897        # define the Iq kernel
898        "#define KERNEL_NAME %s_Iq" % name,
899        call_iq,
900        '#line 1 "%s Iq"' % path,
901        code,
902        clear_iq,
903        "#undef KERNEL_NAME",
904        ]
905
906    iqxy = [
907        # define the Iqxy kernel from the same source with different #defines
908        "#define KERNEL_NAME %s_Iqxy" % name,
909        call_iqxy,
910        '#line 1 "%s Iqxy"' % path,
911        code,
912        clear_iqxy,
913        "#undef KERNEL_NAME",
914    ]
915
916    imagnetic = [
917        # define the Imagnetic kernel
918        "#define KERNEL_NAME %s_Imagnetic" % name,
919        "#define MAGNETIC 1",
920        call_iqxy,
921        '#line 1 "%s Imagnetic"' % path,
922        code,
923        clear_iqxy,
924        "#undef MAGNETIC",
925        "#undef KERNEL_NAME",
926    ]
927
928    return iq, iqxy, imagnetic
929
930
931def load_kernel_module(model_name):
932    # type: (str) -> module
933    """
934    Return the kernel module named in *model_name*.
935
936    If the name ends in *.py* then load it as a custom model using
937    :func:`sasmodels.custom.load_custom_kernel_module`, otherwise
938    load it from :mod:`sasmodels.models`.
939    """
940    if model_name.endswith('.py'):
941        kernel_module = load_custom_kernel_module(model_name)
942    else:
943        try:
944            from sasmodels import models
945            __import__('sasmodels.models.'+model_name)
946            kernel_module = getattr(models, model_name, None)
947        except ImportError:
948            # If the model isn't a built in model, try the plugin directory
949            plugin_path = environ.get('SAS_MODELPATH', None)
950            if plugin_path is not None:
951                file_name = model_name.split(sep)[-1]
952                model_name = plugin_path + sep + file_name + ".py"
953                kernel_module = load_custom_kernel_module(model_name)
954            else:
955                raise
956    return kernel_module
957
958
959section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z'
960                            % re.escape(string.punctuation))
961def _convert_section_titles_to_boldface(lines):
962    # type: (Sequence[str]) -> Iterator[str]
963    """
964    Do the actual work of identifying and converting section headings.
965    """
966    prior = None
967    for line in lines:
968        if prior is None:
969            prior = line
970        elif section_marker.match(line):
971            if len(line) >= len(prior):
972                yield "".join(("**", prior, "**"))
973                prior = None
974            else:
975                yield prior
976                prior = line
977        else:
978            yield prior
979            prior = line
980    if prior is not None:
981        yield prior
982
983
984def convert_section_titles_to_boldface(s):
985    # type: (str) -> str
986    """
987    Use explicit bold-face rather than section headings so that the table of
988    contents is not polluted with section names from the model documentation.
989
990    Sections are identified as the title line followed by a line of punctuation
991    at least as long as the title line.
992    """
993    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
994
995
996def make_doc(model_info):
997    # type: (ModelInfo) -> str
998    """
999    Return the documentation for the model.
1000    """
1001    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
1002    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
1003    docs = model_info.docs if model_info.docs is not None else ""
1004    docs = convert_section_titles_to_boldface(docs)
1005    if model_info.structure_factor:
1006        pars = model_info.parameters.kernel_parameters
1007    else:
1008        pars = (model_info.parameters.common_parameters
1009                + model_info.parameters.kernel_parameters)
1010    partable = make_partable(pars)
1011    subst = dict(id=model_info.id.replace('_', '-'),
1012                 name=model_info.name,
1013                 title=model_info.title,
1014                 parameters=partable,
1015                 returns=Sq_units if model_info.structure_factor else Iq_units,
1016                 docs=docs)
1017    return DOC_HEADER % subst
1018
1019
1020# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog
1021RST_PROLOG = r"""\
1022.. |Ang| unicode:: U+212B
1023.. |Ang^-1| replace:: |Ang|\ :sup:`-1`
1024.. |Ang^2| replace:: |Ang|\ :sup:`2`
1025.. |Ang^-2| replace:: |Ang|\ :sup:`-2`
1026.. |1e-6Ang^-2| replace:: 10\ :sup:`-6`\ |Ang|\ :sup:`-2`
1027.. |Ang^3| replace:: |Ang|\ :sup:`3`
1028.. |Ang^-3| replace:: |Ang|\ :sup:`-3`
1029.. |Ang^-4| replace:: |Ang|\ :sup:`-4`
1030.. |cm^-1| replace:: cm\ :sup:`-1`
1031.. |cm^2| replace:: cm\ :sup:`2`
1032.. |cm^-2| replace:: cm\ :sup:`-2`
1033.. |cm^3| replace:: cm\ :sup:`3`
1034.. |1e15cm^3| replace:: 10\ :sup:`15`\ cm\ :sup:`3`
1035.. |cm^-3| replace:: cm\ :sup:`-3`
1036.. |sr^-1| replace:: sr\ :sup:`-1`
1037
1038.. |cdot| unicode:: U+00B7
1039.. |deg| unicode:: U+00B0
1040.. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3`
1041.. |mg/m^2| replace:: mg\ |cdot|\ m\ :sup:`-2`
1042.. |fm^2| replace:: fm\ :sup:`2`
1043.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1`
1044"""
1045
1046# TODO: make a better fake reference role
1047RST_ROLES = """\
1048.. role:: ref
1049
1050.. role:: numref
1051
1052"""
1053
1054def make_html(model_info):
1055    # type: (ModelInfo) -> str
1056    """
1057    Convert model docs directly to html.
1058    """
1059    from . import rst2html
1060
1061    rst = make_doc(model_info)
1062    return rst2html.rst2html("".join((RST_ROLES, RST_PROLOG, rst)))
1063
1064def view_html(model_name):
1065    # type: (str) -> None
1066    """
1067    Load the model definition and view its help.
1068    """
1069    from . import modelinfo
1070    kernel_module = load_kernel_module(model_name)
1071    info = modelinfo.make_model_info(kernel_module)
1072    view_html_from_info(info)
1073
1074def view_html_from_info(info):
1075    # type: (ModelInfo) -> None
1076    """
1077    View the help for a loaded model definition.
1078    """
1079    from . import rst2html
1080    url = "file://"+dirname(info.filename)+"/"
1081    rst2html.view_html(make_html(info), url=url)
1082
1083def demo_time():
1084    # type: () -> None
1085    """
1086    Show how long it takes to process a model.
1087    """
1088    import datetime
1089    from .modelinfo import make_model_info
1090    from .models import cylinder
1091
1092    tic = datetime.datetime.now()
1093    make_source(make_model_info(cylinder))
1094    toc = (datetime.datetime.now() - tic).total_seconds()
1095    print("time: %g"%toc)
1096
1097
1098def main():
1099    # type: () -> None
1100    """
1101    Program which prints the source produced by the model.
1102    """
1103    from .modelinfo import make_model_info
1104
1105    if len(sys.argv) <= 1:
1106        print("usage: python -m sasmodels.generate modelname")
1107    else:
1108        name = sys.argv[1]
1109        kernel_module = load_kernel_module(name)
1110        model_info = make_model_info(kernel_module)
1111        source = make_source(model_info)
1112        print(source['dll'])
1113
1114
1115if __name__ == "__main__":
1116    main()
Note: See TracBrowser for help on using the repository browser.