source: sasmodels/sasmodels/generate.py @ c036ddb

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

refactor so Iq is not needed if Fq is defined

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