source: sasmodels/sasmodels/generate.py @ 39a06c9

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

Remove references to ER and VR from sasmodels. Refs #1202.

  • Property mode set to 100644
File size: 39.6 KB
Line 
1"""
2SAS model constructor.
3
4Small angle scattering models are defined by a set of kernel functions:
5
6    *Iq(q, p1, p2, ...)* returns the scattering at q for a form with
7    particular dimensions averaged over all orientations.
8
9    *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        m = _FQ_PATTERN.search(code)
706        if m is not None:
707            return True
708    return False
709
710_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE)
711def contains_shell_volume(source):
712    # type: (List[str]) -> bool
713    """
714    Return True if C source defines "void Fq(".
715    """
716    for code in source:
717        m = _SHELL_VOLUME_PATTERN.search(code)
718        if m is not None:
719            return True
720    return False
721
722def _add_source(source, code, path, lineno=1):
723    """
724    Add a file to the list of source code chunks, tagged with path and line.
725    """
726    path = path.replace('\\', '\\\\')
727    source.append('#line %d "%s"' % (lineno, path))
728    source.append(code)
729
730def make_source(model_info):
731    # type: (ModelInfo) -> Dict[str, str]
732    """
733    Generate the OpenCL/ctypes kernel from the module info.
734
735    Uses source files found in the given search path.  Returns None if this
736    is a pure python model, with no C source components.
737    """
738    if callable(model_info.Iq):
739        raise ValueError("can't compile python model")
740        #return None
741
742    # TODO: need something other than volume to indicate dispersion parameters
743    # No volume normalization despite having a volume parameter.
744    # Thickness is labelled a volume in order to trigger polydispersity.
745    # May want a separate dispersion flag, or perhaps a separate category for
746    # disperse, but not volume.  Volume parameters also use relative values
747    # for the distribution rather than the absolute values used by angular
748    # dispersion.  Need to be careful that necessary parameters are available
749    # for computing volume even if we allow non-disperse volume parameters.
750    partable = model_info.parameters
751
752    # Load templates and user code
753    kernel_header = load_template('kernel_header.c')
754    kernel_code = load_template('kernel_iq.c')
755    user_code = [(f, open(f).read()) for f in model_sources(model_info)]
756
757    # Build initial sources
758    source = []
759    _add_source(source, *kernel_header)
760    for path, code in user_code:
761        _add_source(source, code, path)
762    if model_info.c_code:
763        _add_source(source, model_info.c_code, model_info.filename,
764                    lineno=model_info.lineno.get('c_code', 1))
765
766    # Make parameters for q, qx, qy so that we can use them in declarations
767    q, qx, qy, qab, qa, qb, qc \
768        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
769
770    # Generate form_volume function, etc. from body only
771    if isinstance(model_info.form_volume, str):
772        pars = partable.form_volume_parameters
773        source.append(_gen_fn(model_info, 'form_volume', pars))
774    if isinstance(model_info.shell_volume, str):
775        pars = partable.form_volume_parameters
776        source.append(_gen_fn(model_info, 'shell_volume', pars))
777    if isinstance(model_info.Iq, str):
778        pars = [q] + partable.iq_parameters
779        source.append(_gen_fn(model_info, 'Iq', pars))
780    if isinstance(model_info.Iqxy, str):
781        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters
782        source.append(_gen_fn(model_info, 'Iqxy', pars))
783    if isinstance(model_info.Iqac, str):
784        pars = [qab, qc] + partable.iq_parameters
785        source.append(_gen_fn(model_info, 'Iqac', pars))
786    if isinstance(model_info.Iqabc, str):
787        pars = [qa, qb, qc] + partable.iq_parameters
788        source.append(_gen_fn(model_info, 'Iqabc', pars))
789
790    # Check for shell_volume in source
791    is_hollow = contains_shell_volume(source)
792
793    # What kind of 2D model do we need?  Is it consistent with the parameters?
794    xy_mode = find_xy_mode(source)
795    if xy_mode == 'qabc' and not partable.is_asymmetric:
796        raise ValueError("asymmetric oriented models need to define Iqabc")
797    elif xy_mode == 'qac' and partable.is_asymmetric:
798        raise ValueError("symmetric oriented models need to define Iqac")
799    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'):
800        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode)
801    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'):
802        if xy_mode == 'qxy':
803            logger.warn("oriented shapes should define Iqac or Iqabc")
804        else:
805            raise ValueError("Expected function Iqac or Iqabc for oriented shape")
806
807    # Define the parameter table
808    lineno = getframeinfo(currentframe()).lineno + 2
809    source.append('#line %d "sasmodels/generate.py"'%lineno)
810    #source.append('introduce breakage in generate to test lineno reporting')
811    source.append("#define PARAMETER_TABLE \\")
812    source.append("\\\n".join(p.as_definition()
813                              for p in partable.kernel_parameters))
814    # Define the function calls
815    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0"
816    if partable.form_volume_parameters:
817        refs = _call_pars("_v.", partable.form_volume_parameters)
818        if is_hollow:
819            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2)
820        else:
821            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs))
822        if model_info.effective_radius_type:
823            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs))
824    else:
825        # Model doesn't have volume.  We could make the kernel run a little
826        # faster by not using/transferring the volume normalizations, but
827        # the ifdef's reduce readability more than is worthwhile.
828        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)"
829    source.append(call_volume)
830    source.append(call_effective_radius)
831    model_refs = _call_pars("_v.", partable.iq_parameters)
832
833    if model_info.have_Fq:
834        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
835        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars
836        clear_iq = "#undef CALL_FQ"
837    else:
838        pars = ",".join(["_q"] + model_refs)
839        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
840        clear_iq = "#undef CALL_IQ"
841    if xy_mode == 'qabc':
842        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
843        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars
844        clear_iqxy = "#undef CALL_IQ_ABC"
845    elif xy_mode == 'qac':
846        pars = ",".join(["_qa", "_qc"] + model_refs)
847        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
848        clear_iqxy = "#undef CALL_IQ_AC"
849    elif xy_mode == 'qa' and not model_info.have_Fq:
850        pars = ",".join(["_qa"] + model_refs)
851        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
852        clear_iqxy = "#undef CALL_IQ_A"
853    elif xy_mode == 'qa' and model_info.have_Fq:
854        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs)
855        # Note: uses rare C construction (expr1, expr2) which computes
856        # expr1 then expr2 and evaluates to expr2.  This allows us to
857        # leave it looking like a function even though it is returning
858        # its values by reference.
859        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars
860        clear_iqxy = "#undef CALL_FQ_A"
861    elif xy_mode == 'qxy':
862        orientation_refs = _call_pars("_v.", partable.orientation_parameters)
863        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs)
864        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars
865        clear_iqxy = "#undef CALL_IQ_XY"
866        if partable.orientation_parameters:
867            call_iqxy += "\n#define HAVE_THETA"
868            clear_iqxy += "\n#undef HAVE_THETA"
869        if partable.is_asymmetric:
870            call_iqxy += "\n#define HAVE_PSI"
871            clear_iqxy += "\n#undef HAVE_PSI"
872
873
874    magpars = [k-2 for k, p in enumerate(partable.call_parameters)
875               if p.type == 'sld']
876    # Fill in definitions for numbers of parameters
877    source.append("#define MAX_PD %s"%partable.max_pd)
878    source.append("#define NUM_PARS %d"%partable.npars)
879    source.append("#define NUM_VALUES %d" % partable.nvalues)
880    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
881    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
882    source.append("#define PROJECTION %d"%PROJECTION)
883    # TODO: allow mixed python/opencl kernels?
884    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
885    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
886
887    result = {
888        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
889        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
890    }
891    return result
892
893
894def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name):
895    # type: ([str,str], str, str, str) -> List[str]
896    code = kernel[0]
897    path = kernel[1].replace('\\', '\\\\')
898    iq = [
899        # define the Iq kernel
900        "#define KERNEL_NAME %s_Iq" % name,
901        call_iq,
902        '#line 1 "%s Iq"' % path,
903        code,
904        clear_iq,
905        "#undef KERNEL_NAME",
906        ]
907
908    iqxy = [
909        # define the Iqxy kernel from the same source with different #defines
910        "#define KERNEL_NAME %s_Iqxy" % name,
911        call_iqxy,
912        '#line 1 "%s Iqxy"' % path,
913        code,
914        clear_iqxy,
915        "#undef KERNEL_NAME",
916    ]
917
918    imagnetic = [
919        # define the Imagnetic kernel
920        "#define KERNEL_NAME %s_Imagnetic" % name,
921        "#define MAGNETIC 1",
922        call_iqxy,
923        '#line 1 "%s Imagnetic"' % path,
924        code,
925        clear_iqxy,
926        "#undef MAGNETIC",
927        "#undef KERNEL_NAME",
928    ]
929
930    return iq, iqxy, imagnetic
931
932
933def load_kernel_module(model_name):
934    # type: (str) -> module
935    """
936    Return the kernel module named in *model_name*.
937
938    If the name ends in *.py* then load it as a custom model using
939    :func:`sasmodels.custom.load_custom_kernel_module`, otherwise
940    load it from :mod:`sasmodels.models`.
941    """
942    if model_name.endswith('.py'):
943        kernel_module = load_custom_kernel_module(model_name)
944    else:
945        try:
946            from sasmodels import models
947            __import__('sasmodels.models.'+model_name)
948            kernel_module = getattr(models, model_name, None)
949        except ImportError:
950            # If the model isn't a built in model, try the plugin directory
951            plugin_path = environ.get('SAS_MODELPATH', None)
952            if plugin_path is not None:
953                file_name = model_name.split(sep)[-1]
954                model_name = plugin_path + sep + file_name + ".py"
955                kernel_module = load_custom_kernel_module(model_name)
956            else:
957                raise
958    return kernel_module
959
960
961section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z'
962                            % re.escape(string.punctuation))
963def _convert_section_titles_to_boldface(lines):
964    # type: (Sequence[str]) -> Iterator[str]
965    """
966    Do the actual work of identifying and converting section headings.
967    """
968    prior = None
969    for line in lines:
970        if prior is None:
971            prior = line
972        elif section_marker.match(line):
973            if len(line) >= len(prior):
974                yield "".join(("**", prior, "**"))
975                prior = None
976            else:
977                yield prior
978                prior = line
979        else:
980            yield prior
981            prior = line
982    if prior is not None:
983        yield prior
984
985
986def convert_section_titles_to_boldface(s):
987    # type: (str) -> str
988    """
989    Use explicit bold-face rather than section headings so that the table of
990    contents is not polluted with section names from the model documentation.
991
992    Sections are identified as the title line followed by a line of punctuation
993    at least as long as the title line.
994    """
995    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
996
997
998def make_doc(model_info):
999    # type: (ModelInfo) -> str
1000    """
1001    Return the documentation for the model.
1002    """
1003    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
1004    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
1005    docs = model_info.docs if model_info.docs is not None else ""
1006    docs = convert_section_titles_to_boldface(docs)
1007    if model_info.structure_factor:
1008        pars = model_info.parameters.kernel_parameters
1009    else:
1010        pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters
1011    partable = make_partable(pars)
1012    subst = dict(id=model_info.id.replace('_', '-'),
1013                 name=model_info.name,
1014                 title=model_info.title,
1015                 parameters=partable,
1016                 returns=Sq_units if model_info.structure_factor else Iq_units,
1017                 docs=docs)
1018    return DOC_HEADER % subst
1019
1020
1021# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog
1022RST_PROLOG = r"""\
1023.. |Ang| unicode:: U+212B
1024.. |Ang^-1| replace:: |Ang|\ :sup:`-1`
1025.. |Ang^2| replace:: |Ang|\ :sup:`2`
1026.. |Ang^-2| replace:: |Ang|\ :sup:`-2`
1027.. |1e-6Ang^-2| replace:: 10\ :sup:`-6`\ |Ang|\ :sup:`-2`
1028.. |Ang^3| replace:: |Ang|\ :sup:`3`
1029.. |Ang^-3| replace:: |Ang|\ :sup:`-3`
1030.. |Ang^-4| replace:: |Ang|\ :sup:`-4`
1031.. |cm^-1| replace:: cm\ :sup:`-1`
1032.. |cm^2| replace:: cm\ :sup:`2`
1033.. |cm^-2| replace:: cm\ :sup:`-2`
1034.. |cm^3| replace:: cm\ :sup:`3`
1035.. |1e15cm^3| replace:: 10\ :sup:`15`\ cm\ :sup:`3`
1036.. |cm^-3| replace:: cm\ :sup:`-3`
1037.. |sr^-1| replace:: sr\ :sup:`-1`
1038
1039.. |cdot| unicode:: U+00B7
1040.. |deg| unicode:: U+00B0
1041.. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3`
1042.. |mg/m^2| replace:: mg\ |cdot|\ m\ :sup:`-2`
1043.. |fm^2| replace:: fm\ :sup:`2`
1044.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1`
1045"""
1046
1047# TODO: make a better fake reference role
1048RST_ROLES = """\
1049.. role:: ref
1050
1051.. role:: numref
1052
1053"""
1054
1055def make_html(model_info):
1056    # type: (ModelInfo) -> str
1057    """
1058    Convert model docs directly to html.
1059    """
1060    from . import rst2html
1061
1062    rst = make_doc(model_info)
1063    return rst2html.rst2html("".join((RST_ROLES, RST_PROLOG, rst)))
1064
1065def view_html(model_name):
1066    # type: (str) -> None
1067    """
1068    Load the model definition and view its help.
1069    """
1070    from . import modelinfo
1071    kernel_module = load_kernel_module(model_name)
1072    info = modelinfo.make_model_info(kernel_module)
1073    view_html_from_info(info)
1074
1075def view_html_from_info(info):
1076    # type: (ModelInfo) -> None
1077    """
1078    View the help for a loaded model definition.
1079    """
1080    from . import rst2html
1081    url = "file://"+dirname(info.filename)+"/"
1082    rst2html.view_html(make_html(info), url=url)
1083
1084def demo_time():
1085    # type: () -> None
1086    """
1087    Show how long it takes to process a model.
1088    """
1089    import datetime
1090    from .modelinfo import make_model_info
1091    from .models import cylinder
1092
1093    tic = datetime.datetime.now()
1094    make_source(make_model_info(cylinder))
1095    toc = (datetime.datetime.now() - tic).total_seconds()
1096    print("time: %g"%toc)
1097
1098
1099def main():
1100    # type: () -> None
1101    """
1102    Program which prints the source produced by the model.
1103    """
1104    from .modelinfo import make_model_info
1105
1106    if len(sys.argv) <= 1:
1107        print("usage: python -m sasmodels.generate modelname")
1108    else:
1109        name = sys.argv[1]
1110        kernel_module = load_kernel_module(name)
1111        model_info = make_model_info(kernel_module)
1112        source = make_source(model_info)
1113        print(source['dll'])
1114
1115
1116if __name__ == "__main__":
1117    main()
Note: See TracBrowser for help on using the repository browser.