source: sasmodels/sasmodels/generate.py @ b297ba9

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b297ba9 was b297ba9, checked in by Paul Kienzle <pkienzle@…>, 6 months ago

lint

  • Property mode set to 100644
File size: 40.2 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    """
193    Search for the target file relative in the installed application.
194
195    Search first in the location of the generate module in case we are
196    running directly from the distribution.  Search next to the python
197    executable for windows installs.  Search in the ../Resources directory
198    next to the executable for Mac OS/X installs.
199    """
200    path = abspath(dirname(__file__))
201    if exists(joinpath(path, target_file)):
202        return path
203
204    # check next to exe/zip file
205    exepath = dirname(sys.executable)
206    path = joinpath(exepath, external_dir)
207    if exists(joinpath(path, target_file)):
208        return path
209
210    # check in py2app Contents/Resources
211    path = joinpath(exepath, '..', 'Resources', external_dir)
212    if exists(joinpath(path, target_file)):
213        return abspath(path)
214
215    raise RuntimeError('Could not find '+joinpath(external_dir, target_file))
216
217EXTERNAL_DIR = 'sasmodels-data'
218DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_iq.c')
219MODEL_PATH = joinpath(DATA_PATH, 'models')
220
221F16 = np.dtype('float16')
222F32 = np.dtype('float32')
223F64 = np.dtype('float64')
224try:  # CRUFT: older numpy does not support float128
225    F128 = np.dtype('float128')
226except TypeError:
227    F128 = None
228
229# Conversion from units defined in the parameter table for each model
230# to units displayed in the sphinx documentation.
231# This section associates the unit with the macro to use to produce the LaTex
232# code.  The macro itself needs to be defined in sasmodels/doc/rst_prolog.
233#
234# NOTE: there is an RST_PROLOG at the end of this file which is NOT
235# used for the bundled documentation. Still as long as we are defining the
236# macros in two places any new addition should define the macro in both places.
237RST_UNITS = {
238    "Ang": "|Ang|",
239    "1/Ang": "|Ang^-1|",
240    "1/Ang^2": "|Ang^-2|",
241    "Ang^3": "|Ang^3|",
242    "Ang^2": "|Ang^2|",
243    "1e15/cm^3": "|1e15cm^3|",
244    "Ang^3/mol": "|Ang^3|/mol",
245    "1e-6/Ang^2": "|1e-6Ang^-2|",
246    "degrees": "degree",
247    "1/cm": "|cm^-1|",
248    "Ang/cm": "|Ang*cm^-1|",
249    "g/cm^3": "|g/cm^3|",
250    "mg/m^2": "|mg/m^2|",
251    "": "None",
252    }
253
254# Headers for the parameters tables in th sphinx documentation
255PARTABLE_HEADERS = [
256    "Parameter",
257    "Description",
258    "Units",
259    "Default value",
260    ]
261
262# Minimum width for a default value (this is shorter than the column header
263# width, so will be ignored).
264PARTABLE_VALUE_WIDTH = 10
265
266# Documentation header for the module, giving the model name, its short
267# description and its parameter table.  The remainder of the doc comes
268# from the module docstring.
269DOC_HEADER = """.. _%(id)s:
270
271%(name)s
272=======================================================
273
274%(title)s
275
276%(parameters)s
277
278%(returns)s
279
280%(docs)s
281"""
282
283
284def set_integration_size(info, n):
285    # type: (ModelInfo, int) -> None
286    """
287    Update the model definition, replacing the gaussian integration with
288    a gaussian integration of a different size.
289
290    Note: this really ought to be a method in modelinfo, but that leads to
291    import loops.
292    """
293    if info.source and any(lib.startswith('lib/gauss') for lib in info.source):
294        from .gengauss import gengauss
295        path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n)
296        if not exists(path):
297            gengauss(n, path)
298        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')
299                       else lib for lib in info.source]
300
301def format_units(units):
302    # type: (str) -> str
303    """
304    Convert units into ReStructured Text format.
305    """
306    return "string" if isinstance(units, list) else RST_UNITS.get(units, units)
307
308
309def make_partable(pars):
310    # type: (List[Parameter]) -> str
311    """
312    Generate the parameter table to include in the sphinx documentation.
313    """
314    column_widths = [
315        max(len(p.name) for p in pars),
316        max(len(p.description) for p in pars),
317        max(len(format_units(p.units)) for p in pars),
318        PARTABLE_VALUE_WIDTH,
319        ]
320    column_widths = [max(w, len(h))
321                     for w, h in zip(column_widths, PARTABLE_HEADERS)]
322
323    underbar = " ".join("="*w for w in column_widths)
324    lines = [
325        underbar,
326        " ".join("%-*s" % (w, h)
327                 for w, h in zip(column_widths, PARTABLE_HEADERS)),
328        underbar,
329        ]
330    for p in pars:
331        lines.append(" ".join([
332            "%-*s" % (column_widths[0], p.name),
333            "%-*s" % (column_widths[1], p.description),
334            "%-*s" % (column_widths[2], format_units(p.units)),
335            "%*g" % (column_widths[3], p.default),
336            ]))
337    lines.append(underbar)
338    return "\n".join(lines)
339
340
341def _search(search_path, filename):
342    # type: (List[str], str) -> str
343    """
344    Find *filename* in *search_path*.
345
346    Raises ValueError if file does not exist.
347    """
348    for path in search_path:
349        target = joinpath(path, filename)
350        if exists(target):
351            return target
352    raise ValueError("%r not found in %s" % (filename, search_path))
353
354
355def model_sources(model_info):
356    # type: (ModelInfo) -> List[str]
357    """
358    Return a list of the sources file paths for the module.
359    """
360    search_path = [dirname(model_info.filename), MODEL_PATH]
361    return [_search(search_path, f) for f in model_info.source]
362
363
364def dll_timestamp(model_info):
365    # type: (ModelInfo) -> int
366    """
367    Return a timestamp for the model corresponding to the most recently
368    changed file or dependency.
369    """
370    # TODO: fails DRY; templates appear two places.
371    model_templates = [joinpath(DATA_PATH, filename)
372                       for filename in ('kernel_header.c', 'kernel_iq.c')]
373    source_files = (model_sources(model_info)
374                    + model_templates
375                    + [model_info.filename])
376    # Note: file may not exist when it is a standard model from library.zip
377    times = [getmtime(f) for f in source_files if exists(f)]
378    newest = max(times) if times else 0
379    return newest
380
381def ocl_timestamp(model_info):
382    # type: (ModelInfo) -> int
383    """
384    Return a timestamp for the model corresponding to the most recently
385    changed file or dependency.
386
387    Note that this does not look at the time stamps for the OpenCL header
388    information since that need not trigger a recompile of the DLL.
389    """
390    # TODO: fails DRY; templates appear two places.
391    model_templates = [joinpath(DATA_PATH, filename)
392                       for filename in ('kernel_header.c', 'kernel_iq.c')]
393    source_files = (model_sources(model_info)
394                    + model_templates
395                    + [model_info.filename])
396    # Note: file may not exist when it is a standard model from library.zip
397    times = [getmtime(f) for f in source_files if exists(f)]
398    newest = max(times) if times else 0
399    return newest
400
401def tag_source(source):
402    # type: (str) -> str
403    """
404    Return a unique tag for the source code.
405    """
406    # Note: need 0xffffffff&val to force an unsigned 32-bit number
407    try:
408        source = source.encode('utf8')
409    except AttributeError: # bytes has no encode attribute in python 3
410        pass
411    return "%08X"%(0xffffffff&crc32(source))
412
413def convert_type(source, dtype):
414    # type: (str, np.dtype) -> str
415    """
416    Convert code from double precision to the desired type.
417
418    Floating point constants are tagged with 'f' for single precision or 'L'
419    for long double precision.
420    """
421    source = _fix_tgmath_int(source)
422    if dtype == F16:
423        fbytes = 2
424        source = _convert_type(source, "half", "f")
425    elif dtype == F32:
426        fbytes = 4
427        source = _convert_type(source, "float", "f")
428    elif dtype == F64:
429        fbytes = 8
430        # no need to convert if it is already double
431    elif dtype == F128:
432        fbytes = 16
433        source = _convert_type(source, "long double", "L")
434    else:
435        raise ValueError("Unexpected dtype in source conversion: %s" % dtype)
436    return ("#define FLOAT_SIZE %d\n" % fbytes)+source
437
438
439def _convert_type(source, type_name, constant_flag):
440    # type: (str, str, str) -> str
441    """
442    Replace 'double' with *type_name* in *source*, tagging floating point
443    constants with *constant_flag*.
444    """
445    # Convert double keyword to float/long double/half.
446    # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16.
447    # Assume complex numbers are represented as cdouble which is typedef'd
448    # to double2.
449    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
450                    r'\1%s\2'%type_name, source)
451    source = _tag_float(source, constant_flag)
452    return source
453
454TGMATH_INT_RE = re.compile(r"""
455(?: # Non-capturing match; not lookbehind since pattern length is variable
456  \b              # word boundary
457   # various math functions
458  (a?(sin|cos|tan)h? | atan2
459   | erfc? | tgamma
460   | exp(2|10|m1)? | log(2|10|1p)? | pow[nr]? | sqrt | rsqrt | rootn
461   | fabs | fmax | fmin
462   )
463  \s*[(]\s*       # open parenthesis
464)
465[+-]?(0|[1-9]\d*) # integer
466(?=               # lookahead match: don't want to move from end of int
467  \s*[,)]         # comma or close parenthesis for end of argument
468)                 # end lookahead
469""", re.VERBOSE)
470def _fix_tgmath_int(source):
471    # type: (str) -> str
472    """
473    Replace f(integer) with f(integer.) for sin, cos, pow, etc.
474
475    OS X OpenCL complains that it can't resolve the type generic calls to
476    the standard math functions when they are called with integer constants,
477    but this does not happen with the Windows Intel driver for example.
478    To avoid confusion on the matrix marketplace, automatically promote
479    integers to floats if we recognize them in the source.
480
481    The specific functions we look for are:
482
483        trigonometric: sin, asin, sinh, asinh, etc., and atan2
484        exponential:   exp, exp2, exp10, expm1, log, log2, log10, logp1
485        power:         pow, pown, powr, sqrt, rsqrt, rootn
486        special:       erf, erfc, tgamma
487        float:         fabs, fmin, fmax
488
489    Note that we don't convert the second argument of dual argument
490    functions: atan2, fmax, fmin, pow, powr.  This could potentially
491    be a problem for pow(x, 2), but that case seems to work without change.
492    """
493    out = TGMATH_INT_RE.sub(r'\g<0>.', source)
494    return out
495
496
497# Floating point regular expression
498#
499# Define parts:
500#
501#    E = [eE][+-]?\d+    : Exponent
502#    P = [.]             : Decimal separator
503#    N = [1-9]\d*        : Natural number, no leading zeros
504#    Z = 0               : Zero
505#    F = \d+             : Fractional number, maybe leading zeros
506#    F? = \d*            : Optional fractional number
507#
508# We want to reject bare natural numbers and bare decimal points, so we
509# need to tediously outline the cases where we have either a fraction or
510# an exponent:
511#
512#   ( ZP | ZPF | ZE | ZPE | ZPFE | NP | NPF | NE | NPE | NPFE | PF | PFE )
513#
514#
515# We can then join cases by making parts optional.  The following are
516# some ways to do this:
517#
518#   ( (Z|N)(P|PF|E|PE|PFE) | PFE? )                   # Split on lead
519#     => ( (Z|N)(PF?|(PF?)?E) | PFE? )
520#   ( ((Z|N)PF?|PF)E? | (Z|N)E)                       # Split on point
521#   ( (ZP|ZPF|NP|NPF|PF) | (Z|ZP|ZPF|N|NP|NPF|PF)E )  # Split on E
522#     => ( ((Z|N)PF?|PF) | ((Z|N)(PF?)? | PF) E )
523FLOAT_RE = re.compile(r"""
524    (?<!\w)  # use negative lookbehind since '.' confuses \b test
525    # use split on lead to match float ( (Z|N)(PF?|(PF?)?E) | PFE? )
526    ( ( 0 | [1-9]\d* )                     # ( ( Z | N )
527      ([.]\d* | ([.]\d*)? [eE][+-]?\d+ )   #   (PF? | (PF?)? E )
528    | [.]\d+ ([eE][+-]?\d+)?               # | PF (E)?
529    )                                      # )
530    (?!\w)  # use negative lookahead since '.' confuses \b test
531    """, re.VERBOSE)
532def _tag_float(source, constant_flag):
533    # Convert floating point constants to single by adding 'f' to the end,
534    # or long double with an 'L' suffix.  OS/X complains if you don't do this.
535    out = FLOAT_RE.sub(r'\g<0>%s'%constant_flag, source)
536    #print("in",repr(source),"out",repr(out), constant_flag)
537    return out
538
539def test_tag_float():
540    """Check that floating point constants are identified and tagged with 'f'"""
541
542    cases = """
543ZP  : 0.
544ZPF : 0.0,0.01,0.1
545Z  E: 0e+001
546ZP E: 0.E0
547ZPFE: 0.13e-031
548NP  : 1., 12.
549NPF : 1.0001, 1.1, 1.0
550N  E: 1e0, 37E-080
551NP E: 1.e0, 37.E-080
552NPFE: 845.017e+22
553 PF : .1, .0, .0100
554 PFE: .6e+9, .82E-004
555# isolated cases
5560.
5571e0
5580.13e-013
559# untouched
560struct3.e3, 03.05.67, 37
561# expressions
5623.75+-1.6e-7-27+13.2
563a3.e2 - 0.
5644*atan(1)
5654.*atan(1.)
566"""
567
568    output = """
569ZP  : 0.f
570ZPF : 0.0f,0.01f,0.1f
571Z  E: 0e+001f
572ZP E: 0.E0f
573ZPFE: 0.13e-031f
574NP  : 1.f, 12.f
575NPF : 1.0001f, 1.1f, 1.0f
576N  E: 1e0f, 37E-080f
577NP E: 1.e0f, 37.E-080f
578NPFE: 845.017e+22f
579 PF : .1f, .0f, .0100f
580 PFE: .6e+9f, .82E-004f
581# isolated cases
5820.f
5831e0f
5840.13e-013f
585# untouched
586struct3.e3, 03.05.67, 37
587# expressions
5883.75f+-1.6e-7f-27+13.2f
589a3.e2 - 0.f
5904*atan(1)
5914.f*atan(1.f)
592"""
593
594    for case_in, case_out in zip(cases.split('\n'), output.split('\n')):
595        out = _tag_float(case_in, 'f')
596        assert case_out == out, "%r => %r"%(case_in, out)
597
598
599def kernel_name(model_info, variant):
600    # type: (ModelInfo, str) -> str
601    """
602    Name of the exported kernel symbol.
603
604    *variant* is "Iq", "Iqxy" or "Imagnetic".
605    """
606    return model_info.name + "_" + variant
607
608
609def indent(s, depth):
610    # type: (str, int) -> str
611    """
612    Indent a string of text with *depth* additional spaces on each line.
613    """
614    spaces = " "*depth
615    interline_separator = "\n" + spaces
616    return spaces + interline_separator.join(s.split("\n"))
617
618
619_template_cache = {}  # type: Dict[str, Tuple[int, str, str]]
620def load_template(filename):
621    # type: (str) -> str
622    """
623    Load template file from sasmodels resource directory.
624    """
625    path = joinpath(DATA_PATH, filename)
626    mtime = getmtime(path)
627    if filename not in _template_cache or mtime > _template_cache[filename][0]:
628        with open(path) as fid:
629            _template_cache[filename] = (mtime, fid.read(), path)
630    return _template_cache[filename][1], path
631
632
633_FN_TEMPLATE = """\
634double %(name)s(%(pars)s);
635double %(name)s(%(pars)s) {
636#line %(line)d "%(filename)s"
637    %(body)s
638}
639
640"""
641def _gen_fn(model_info, name, pars):
642    # type: (ModelInfo, str, List[Parameter]) -> str
643    """
644    Generate a function given pars and body.
645
646    Returns the following string::
647
648         double fn(double a, double b, ...);
649         double fn(double a, double b, ...) {
650             ....
651         }
652    """
653    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
654    body = getattr(model_info, name)
655    filename = model_info.filename
656    # Note: if symbol is defined strangely in the module then default it to 1
657    lineno = model_info.lineno.get(name, 1)
658    return _FN_TEMPLATE % {
659        'name': name, 'pars': par_decl, 'body': body,
660        'filename': filename.replace('\\', '\\\\'), 'line': lineno,
661    }
662
663
664def _call_pars(prefix, pars):
665    # type: (str, List[Parameter]) -> List[str]
666    """
667    Return a list of *prefix+parameter* from parameter items.
668
669    *prefix* should be "v." if v is a struct.
670    """
671    return [p.as_call_reference(prefix) for p in pars]
672
673
674_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]",
675                           flags=re.MULTILINE)
676def find_xy_mode(source):
677    # type: (List[str]) -> bool
678    """
679    Return the xy mode as qa, qac, qabc or qxy.
680
681    Note this is not a C parser, and so can be easily confused by
682    non-standard syntax.  Also, it will incorrectly identify the following
683    as having 2D models::
684
685        /*
686        double Iqac(qab, qc, ...) { ... fill this in later ... }
687        */
688
689    If you want to comment out the function, use // on the front of the
690    line::
691
692        /*
693        // double Iqac(qab, qc, ...) { ... fill this in later ... }
694        */
695
696    """
697    for code in source:
698        m = _IQXY_PATTERN.search(code)
699        if m is not None:
700            return m.group('mode')
701    return 'qa'
702
703# Note: not presently used.  Need to know whether Fq is available before
704# trying to compile the source.  Leave the code here in case we decide that
705# define have_Fq for each form factor is too tedious and error prone.
706_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE)
707def contains_Fq(source):
708    # type: (List[str]) -> bool
709    """
710    Return True if C source defines "void Fq(".
711    """
712    for code in source:
713        if _FQ_PATTERN.search(code) is not None:
714            return True
715    return False
716
717_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]",
718                                   flags=re.MULTILINE)
719def contains_shell_volume(source):
720    # type: (List[str]) -> bool
721    """
722    Return True if C source defines "double shell_volume(".
723    """
724    for code in source:
725        if _SHELL_VOLUME_PATTERN.search(code) is not None:
726            return True
727    return False
728
729def _add_source(source, code, path, lineno=1):
730    """
731    Add a file to the list of source code chunks, tagged with path and line.
732    """
733    path = path.replace('\\', '\\\\')
734    source.append('#line %d "%s"' % (lineno, path))
735    source.append(code)
736
737def make_source(model_info):
738    # type: (ModelInfo) -> Dict[str, str]
739    """
740    Generate the OpenCL/ctypes kernel from the module info.
741
742    Uses source files found in the given search path.  Returns None if this
743    is a pure python model, with no C source components.
744    """
745    if callable(model_info.Iq):
746        raise ValueError("can't compile python model")
747        #return None
748
749    # TODO: need something other than volume to indicate dispersion parameters
750    # No volume normalization despite having a volume parameter.
751    # Thickness is labelled a volume in order to trigger polydispersity.
752    # May want a separate dispersion flag, or perhaps a separate category for
753    # disperse, but not volume.  Volume parameters also use relative values
754    # for the distribution rather than the absolute values used by angular
755    # dispersion.  Need to be careful that necessary parameters are available
756    # for computing volume even if we allow non-disperse volume parameters.
757    partable = model_info.parameters
758
759    # Load templates and user code
760    kernel_header = load_template('kernel_header.c')
761    kernel_code = load_template('kernel_iq.c')
762    user_code = [(f, open(f).read()) for f in model_sources(model_info)]
763
764    # Build initial sources
765    source = []
766    _add_source(source, *kernel_header)
767    for path, code in user_code:
768        _add_source(source, code, path)
769    if model_info.c_code:
770        _add_source(source, model_info.c_code, model_info.filename,
771                    lineno=model_info.lineno.get('c_code', 1))
772
773    # Make parameters for q, qx, qy so that we can use them in declarations
774    q, qx, qy, qab, qa, qb, qc \
775        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
776
777    # Generate form_volume function, etc. from body only
778    if isinstance(model_info.form_volume, str):
779        pars = partable.form_volume_parameters
780        source.append(_gen_fn(model_info, 'form_volume', pars))
781    if isinstance(model_info.shell_volume, str):
782        pars = partable.form_volume_parameters
783        source.append(_gen_fn(model_info, 'shell_volume', pars))
784    if isinstance(model_info.Iq, str):
785        pars = [q] + partable.iq_parameters
786        source.append(_gen_fn(model_info, 'Iq', pars))
787    if isinstance(model_info.Iqxy, str):
788        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters
789        source.append(_gen_fn(model_info, 'Iqxy', pars))
790    if isinstance(model_info.Iqac, str):
791        pars = [qab, qc] + partable.iq_parameters
792        source.append(_gen_fn(model_info, 'Iqac', pars))
793    if isinstance(model_info.Iqabc, str):
794        pars = [qa, qb, qc] + partable.iq_parameters
795        source.append(_gen_fn(model_info, 'Iqabc', pars))
796
797    # Check for shell_volume in source
798    is_hollow = contains_shell_volume(source)
799
800    # What kind of 2D model do we need?  Is it consistent with the parameters?
801    xy_mode = find_xy_mode(source)
802    if xy_mode == 'qabc' and not partable.is_asymmetric:
803        raise ValueError("asymmetric oriented models need to define Iqabc")
804    elif xy_mode == 'qac' and partable.is_asymmetric:
805        raise ValueError("symmetric oriented models need to define Iqac")
806    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'):
807        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode)
808    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'):
809        if xy_mode == 'qxy':
810            logger.warn("oriented shapes should define Iqac or Iqabc")
811        else:
812            raise ValueError("Expected Iqac or Iqabc for oriented shape")
813
814    # Define the parameter table
815    lineno = getframeinfo(currentframe()).lineno + 2
816    source.append('#line %d "sasmodels/generate.py"'%lineno)
817    #source.append('introduce breakage in generate to test lineno reporting')
818    source.append("#define PARAMETER_TABLE \\")
819    source.append("\\\n".join(p.as_definition()
820                              for p in partable.kernel_parameters))
821    # Define the function calls
822    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0"
823    if partable.form_volume_parameters:
824        refs = _call_pars("_v.", partable.form_volume_parameters)
825        if is_hollow:
826            call_volume = (
827                "#define CALL_VOLUME(_form, _shell, _v) "
828                "do { _form = form_volume(%s); _shell = shell_volume(%s); } "
829                "while (0)") % ((",".join(refs),)*2)
830        else:
831            call_volume = (
832                "#define CALL_VOLUME(_form, _shell, _v) "
833                "do { _form = _shell = form_volume(%s); } "
834                "while (0)") % (",".join(refs))
835        if model_info.effective_radius_type:
836            call_effective_radius = (
837                "#define CALL_EFFECTIVE_RADIUS(_mode, _v) "
838                "effective_radius(_mode, %s)") % (",".join(refs))
839    else:
840        # Model doesn't have volume.  We could make the kernel run a little
841        # faster by not using/transferring the volume normalizations, but
842        # the ifdef's reduce readability more than is worthwhile.
843        call_volume = (
844            "#define CALL_VOLUME(_form, _shell, _v) "
845            "do { _form = _shell = 1.0; } while (0)")
846    source.append(call_volume)
847    source.append(call_effective_radius)
848    model_refs = _call_pars("_v.", partable.iq_parameters)
849
850    if model_info.have_Fq:
851        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
852        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars
853        clear_iq = "#undef CALL_FQ"
854    else:
855        pars = ",".join(["_q"] + model_refs)
856        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
857        clear_iq = "#undef CALL_IQ"
858    if xy_mode == 'qabc':
859        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
860        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars
861        clear_iqxy = "#undef CALL_IQ_ABC"
862    elif xy_mode == 'qac':
863        pars = ",".join(["_qa", "_qc"] + model_refs)
864        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
865        clear_iqxy = "#undef CALL_IQ_AC"
866    elif xy_mode == 'qa' and not model_info.have_Fq:
867        pars = ",".join(["_qa"] + model_refs)
868        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
869        clear_iqxy = "#undef CALL_IQ_A"
870    elif xy_mode == 'qa' and model_info.have_Fq:
871        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs)
872        # Note: uses rare C construction (expr1, expr2) which computes
873        # expr1 then expr2 and evaluates to expr2.  This allows us to
874        # leave it looking like a function even though it is returning
875        # its values by reference.
876        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars
877        clear_iqxy = "#undef CALL_FQ_A"
878    elif xy_mode == 'qxy':
879        orientation_refs = _call_pars("_v.", partable.orientation_parameters)
880        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs)
881        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars
882        clear_iqxy = "#undef CALL_IQ_XY"
883        if partable.orientation_parameters:
884            call_iqxy += "\n#define HAVE_THETA"
885            clear_iqxy += "\n#undef HAVE_THETA"
886        if partable.is_asymmetric:
887            call_iqxy += "\n#define HAVE_PSI"
888            clear_iqxy += "\n#undef HAVE_PSI"
889
890
891    magpars = [k-2 for k, p in enumerate(partable.call_parameters)
892               if p.type == 'sld']
893    # Fill in definitions for numbers of parameters
894    source.append("#define MAX_PD %s"%partable.max_pd)
895    source.append("#define NUM_PARS %d"%partable.npars)
896    source.append("#define NUM_VALUES %d" % partable.nvalues)
897    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
898    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
899    source.append("#define PROJECTION %d"%PROJECTION)
900    wrappers = _kernels(kernel_code, call_iq, clear_iq,
901                        call_iqxy, clear_iqxy, model_info.name)
902    code = '\n'.join(source + wrappers[0] + wrappers[1] + wrappers[2])
903
904    # Note: Identical code for dll and opencl.  This may not always be the case
905    # so leave them as separate entries in the returned value.
906    result = {'dll': code, 'opencl': code}
907    return result
908
909
910def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name):
911    # type: ([str,str], str, str, str) -> List[str]
912    code = kernel[0]
913    path = kernel[1].replace('\\', '\\\\')
914    iq = [
915        # define the Iq kernel
916        "#define KERNEL_NAME %s_Iq" % name,
917        call_iq,
918        '#line 1 "%s Iq"' % path,
919        code,
920        clear_iq,
921        "#undef KERNEL_NAME",
922        ]
923
924    iqxy = [
925        # define the Iqxy kernel from the same source with different #defines
926        "#define KERNEL_NAME %s_Iqxy" % name,
927        call_iqxy,
928        '#line 1 "%s Iqxy"' % path,
929        code,
930        clear_iqxy,
931        "#undef KERNEL_NAME",
932    ]
933
934    imagnetic = [
935        # define the Imagnetic kernel
936        "#define KERNEL_NAME %s_Imagnetic" % name,
937        "#define MAGNETIC 1",
938        call_iqxy,
939        '#line 1 "%s Imagnetic"' % path,
940        code,
941        clear_iqxy,
942        "#undef MAGNETIC",
943        "#undef KERNEL_NAME",
944    ]
945
946    return iq, iqxy, imagnetic
947
948
949def load_kernel_module(model_name):
950    # type: (str) -> module
951    """
952    Return the kernel module named in *model_name*.
953
954    If the name ends in *.py* then load it as a custom model using
955    :func:`sasmodels.custom.load_custom_kernel_module`, otherwise
956    load it from :mod:`sasmodels.models`.
957    """
958    if model_name.endswith('.py'):
959        kernel_module = load_custom_kernel_module(model_name)
960    else:
961        try:
962            from sasmodels import models
963            __import__('sasmodels.models.'+model_name)
964            kernel_module = getattr(models, model_name, None)
965        except ImportError:
966            # If the model isn't a built in model, try the plugin directory
967            plugin_path = environ.get('SAS_MODELPATH', None)
968            if plugin_path is not None:
969                file_name = model_name.split(sep)[-1]
970                model_name = plugin_path + sep + file_name + ".py"
971                kernel_module = load_custom_kernel_module(model_name)
972            else:
973                raise
974    return kernel_module
975
976
977section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z'
978                            % re.escape(string.punctuation))
979def _convert_section_titles_to_boldface(lines):
980    # type: (Sequence[str]) -> Iterator[str]
981    """
982    Do the actual work of identifying and converting section headings.
983    """
984    prior = None
985    for line in lines:
986        if prior is None:
987            prior = line
988        elif section_marker.match(line):
989            if len(line) >= len(prior):
990                yield "".join(("**", prior, "**"))
991                prior = None
992            else:
993                yield prior
994                prior = line
995        else:
996            yield prior
997            prior = line
998    if prior is not None:
999        yield prior
1000
1001
1002def convert_section_titles_to_boldface(s):
1003    # type: (str) -> str
1004    """
1005    Use explicit bold-face rather than section headings so that the table of
1006    contents is not polluted with section names from the model documentation.
1007
1008    Sections are identified as the title line followed by a line of punctuation
1009    at least as long as the title line.
1010    """
1011    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
1012
1013
1014def make_doc(model_info):
1015    # type: (ModelInfo) -> str
1016    """
1017    Return the documentation for the model.
1018    """
1019    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
1020    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
1021    docs = model_info.docs if model_info.docs is not None else ""
1022    docs = convert_section_titles_to_boldface(docs)
1023    if model_info.structure_factor:
1024        pars = model_info.parameters.kernel_parameters
1025    else:
1026        pars = (model_info.parameters.common_parameters
1027                + model_info.parameters.kernel_parameters)
1028    partable = make_partable(pars)
1029    subst = dict(id=model_info.id.replace('_', '-'),
1030                 name=model_info.name,
1031                 title=model_info.title,
1032                 parameters=partable,
1033                 returns=Sq_units if model_info.structure_factor else Iq_units,
1034                 docs=docs)
1035    return DOC_HEADER % subst
1036
1037
1038# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog
1039RST_PROLOG = r"""\
1040.. |Ang| unicode:: U+212B
1041.. |Ang^-1| replace:: |Ang|\ :sup:`-1`
1042.. |Ang^2| replace:: |Ang|\ :sup:`2`
1043.. |Ang^-2| replace:: |Ang|\ :sup:`-2`
1044.. |1e-6Ang^-2| replace:: 10\ :sup:`-6`\ |Ang|\ :sup:`-2`
1045.. |Ang^3| replace:: |Ang|\ :sup:`3`
1046.. |Ang^-3| replace:: |Ang|\ :sup:`-3`
1047.. |Ang^-4| replace:: |Ang|\ :sup:`-4`
1048.. |cm^-1| replace:: cm\ :sup:`-1`
1049.. |cm^2| replace:: cm\ :sup:`2`
1050.. |cm^-2| replace:: cm\ :sup:`-2`
1051.. |cm^3| replace:: cm\ :sup:`3`
1052.. |1e15cm^3| replace:: 10\ :sup:`15`\ cm\ :sup:`3`
1053.. |cm^-3| replace:: cm\ :sup:`-3`
1054.. |sr^-1| replace:: sr\ :sup:`-1`
1055
1056.. |cdot| unicode:: U+00B7
1057.. |deg| unicode:: U+00B0
1058.. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3`
1059.. |mg/m^2| replace:: mg\ |cdot|\ m\ :sup:`-2`
1060.. |fm^2| replace:: fm\ :sup:`2`
1061.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1`
1062"""
1063
1064# TODO: make a better fake reference role
1065RST_ROLES = """\
1066.. role:: ref
1067
1068.. role:: numref
1069
1070"""
1071
1072def make_html(model_info):
1073    # type: (ModelInfo) -> str
1074    """
1075    Convert model docs directly to html.
1076    """
1077    from . import rst2html
1078
1079    rst = make_doc(model_info)
1080    return rst2html.rst2html("".join((RST_ROLES, RST_PROLOG, rst)))
1081
1082def view_html(model_name):
1083    # type: (str) -> None
1084    """
1085    Load the model definition and view its help.
1086    """
1087    from . import modelinfo
1088    kernel_module = load_kernel_module(model_name)
1089    info = modelinfo.make_model_info(kernel_module)
1090    view_html_from_info(info)
1091
1092def view_html_from_info(info):
1093    # type: (ModelInfo) -> None
1094    """
1095    View the help for a loaded model definition.
1096    """
1097    from . import rst2html
1098    url = "file://"+dirname(info.filename)+"/"
1099    rst2html.view_html(make_html(info), url=url)
1100
1101def demo_time():
1102    # type: () -> None
1103    """
1104    Show how long it takes to process a model.
1105    """
1106    import datetime
1107    from .modelinfo import make_model_info
1108    from .models import cylinder
1109
1110    tic = datetime.datetime.now()
1111    make_source(make_model_info(cylinder))
1112    toc = (datetime.datetime.now() - tic).total_seconds()
1113    print("time: %g"%toc)
1114
1115
1116def main():
1117    # type: () -> None
1118    """
1119    Program which prints the source produced by the model.
1120    """
1121    from .modelinfo import make_model_info
1122
1123    if len(sys.argv) <= 1:
1124        print("usage: python -m sasmodels.generate modelname")
1125    else:
1126        name = sys.argv[1]
1127        kernel_module = load_kernel_module(name)
1128        model_info = make_model_info(kernel_module)
1129        source = make_source(model_info)
1130        print(source['dll'])
1131
1132
1133if __name__ == "__main__":
1134    main()
Note: See TracBrowser for help on using the repository browser.