source: sasmodels/sasmodels/generate.py @ e44432d

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

support hollow models in structure factor calculations

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