source: sasmodels/sasmodels/generate.py @ 6cbdcd4

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

code cleanup: kernel_template.c/kernel_iq.cl are no longer used by generate

  • Property mode set to 100644
File size: 36.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.path import abspath, dirname, join as joinpath, exists, getmtime
172import re
173import string
174from zlib import crc32
175from inspect import currentframe, getframeinfo
176import logging
177
178import numpy as np  # type: ignore
179
180from .modelinfo import Parameter
181from .custom import load_custom_kernel_module
182
183# pylint: disable=unused-import
184try:
185    from typing import Tuple, Sequence, Iterator, Dict
186    from .modelinfo import ModelInfo
187except ImportError:
188    pass
189# pylint: enable=unused-import
190
191logger = logging.getLogger(__name__)
192
193# jitter projection to use in the kernel code.  See explore/jitter.py
194# for details.  To change it from a program, set generate.PROJECTION.
195PROJECTION = 1
196
197def get_data_path(external_dir, target_file):
198    path = abspath(dirname(__file__))
199    if exists(joinpath(path, target_file)):
200        return path
201
202    # check next to exe/zip file
203    exepath = dirname(sys.executable)
204    path = joinpath(exepath, external_dir)
205    if exists(joinpath(path, target_file)):
206        return path
207
208    # check in py2app Contents/Resources
209    path = joinpath(exepath, '..', 'Resources', external_dir)
210    if exists(joinpath(path, target_file)):
211        return abspath(path)
212
213    raise RuntimeError('Could not find '+joinpath(external_dir, target_file))
214
215EXTERNAL_DIR = 'sasmodels-data'
216DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_iq.c')
217MODEL_PATH = joinpath(DATA_PATH, 'models')
218
219F16 = np.dtype('float16')
220F32 = np.dtype('float32')
221F64 = np.dtype('float64')
222try:  # CRUFT: older numpy does not support float128
223    F128 = np.dtype('float128')
224except TypeError:
225    F128 = None
226
227# Conversion from units defined in the parameter table for each model
228# to units displayed in the sphinx documentation.
229# This section associates the unit with the macro to use to produce the LaTex
230# code.  The macro itself needs to be defined in sasmodels/doc/rst_prolog.
231#
232# NOTE: there is an RST_PROLOG at the end of this file which is NOT
233# used for the bundled documentation. Still as long as we are defining the macros
234# in two places any new addition should define the macro in both places.
235RST_UNITS = {
236    "Ang": "|Ang|",
237    "1/Ang": "|Ang^-1|",
238    "1/Ang^2": "|Ang^-2|",
239    "Ang^3": "|Ang^3|",
240    "Ang^2": "|Ang^2|",
241    "1e15/cm^3": "|1e15cm^3|",
242    "Ang^3/mol": "|Ang^3|/mol",
243    "1e-6/Ang^2": "|1e-6Ang^-2|",
244    "degrees": "degree",
245    "1/cm": "|cm^-1|",
246    "Ang/cm": "|Ang*cm^-1|",
247    "g/cm^3": "|g/cm^3|",
248    "mg/m^2": "|mg/m^2|",
249    "": "None",
250    }
251
252# Headers for the parameters tables in th sphinx documentation
253PARTABLE_HEADERS = [
254    "Parameter",
255    "Description",
256    "Units",
257    "Default value",
258    ]
259
260# Minimum width for a default value (this is shorter than the column header
261# width, so will be ignored).
262PARTABLE_VALUE_WIDTH = 10
263
264# Documentation header for the module, giving the model name, its short
265# description and its parameter table.  The remainder of the doc comes
266# from the module docstring.
267DOC_HEADER = """.. _%(id)s:
268
269%(name)s
270=======================================================
271
272%(title)s
273
274%(parameters)s
275
276%(returns)s
277
278%(docs)s
279"""
280
281
282def set_integration_size(info, n):
283    # type: (ModelInfo, int) -> None
284    """
285    Update the model definition, replacing the gaussian integration with
286    a gaussian integration of a different size.
287
288    Note: this really ought to be a method in modelinfo, but that leads to
289    import loops.
290    """
291    if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)):
292        import os.path
293        from .gengauss import gengauss
294        path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n)
295        if not os.path.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    sep = " ".join("="*w for w in column_widths)
323    lines = [
324        sep,
325        " ".join("%-*s" % (w, h)
326                 for w, h in zip(column_widths, PARTABLE_HEADERS)),
327        sep,
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(sep)
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    sep = "\n" + spaces
615    return spaces + sep.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    path = joinpath(DATA_PATH, filename)
622    mtime = getmtime(path)
623    if filename not in _template_cache or mtime > _template_cache[filename][0]:
624        with open(path) as fid:
625            _template_cache[filename] = (mtime, fid.read(), path)
626    return _template_cache[filename][1], path
627
628
629_FN_TEMPLATE = """\
630double %(name)s(%(pars)s);
631double %(name)s(%(pars)s) {
632#line %(line)d "%(filename)s"
633    %(body)s
634}
635
636"""
637def _gen_fn(model_info, name, pars):
638    # type: (ModelInfo, str, List[Parameter]) -> str
639    """
640    Generate a function given pars and body.
641
642    Returns the following string::
643
644         double fn(double a, double b, ...);
645         double fn(double a, double b, ...) {
646             ....
647         }
648    """
649    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
650    body = getattr(model_info, name)
651    filename = model_info.filename
652    # Note: if symbol is defined strangely in the module then default it to 1
653    lineno = model_info.lineno.get(name, 1)
654    return _FN_TEMPLATE % {
655        'name': name, 'pars': par_decl, 'body': body,
656        'filename': filename.replace('\\', '\\\\'), 'line': lineno,
657    }
658
659
660def _call_pars(prefix, pars):
661    # type: (str, List[Parameter]) -> List[str]
662    """
663    Return a list of *prefix+parameter* from parameter items.
664
665    *prefix* should be "v." if v is a struct.
666    """
667    return [p.as_call_reference(prefix) for p in pars]
668
669
670# type in IQXY pattern could be single, float, double, long double, ...
671_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]",
672                           flags=re.MULTILINE)
673def find_xy_mode(source):
674    # type: (List[str]) -> bool
675    """
676    Return the xy mode as qa, qac, qabc or qxy.
677
678    Note this is not a C parser, and so can be easily confused by
679    non-standard syntax.  Also, it will incorrectly identify the following
680    as having 2D models::
681
682        /*
683        double Iqac(qab, qc, ...) { ... fill this in later ... }
684        */
685
686    If you want to comment out the function, use // on the front of the
687    line::
688
689        /*
690        // double Iqac(qab, qc, ...) { ... fill this in later ... }
691        */
692
693    """
694    for code in source:
695        m = _IQXY_PATTERN.search(code)
696        if m is not None:
697            return m.group('mode')
698    return 'qa'
699
700
701def _add_source(source, code, path, lineno=1):
702    """
703    Add a file to the list of source code chunks, tagged with path and line.
704    """
705    path = path.replace('\\', '\\\\')
706    source.append('#line %d "%s"' % (lineno, path))
707    source.append(code)
708
709def make_source(model_info):
710    # type: (ModelInfo) -> Dict[str, str]
711    """
712    Generate the OpenCL/ctypes kernel from the module info.
713
714    Uses source files found in the given search path.  Returns None if this
715    is a pure python model, with no C source components.
716    """
717    if callable(model_info.Iq):
718        raise ValueError("can't compile python model")
719        #return None
720
721    # TODO: need something other than volume to indicate dispersion parameters
722    # No volume normalization despite having a volume parameter.
723    # Thickness is labelled a volume in order to trigger polydispersity.
724    # May want a separate dispersion flag, or perhaps a separate category for
725    # disperse, but not volume.  Volume parameters also use relative values
726    # for the distribution rather than the absolute values used by angular
727    # dispersion.  Need to be careful that necessary parameters are available
728    # for computing volume even if we allow non-disperse volume parameters.
729
730    partable = model_info.parameters
731
732    # Load templates and user code
733    kernel_header = load_template('kernel_header.c')
734    kernel_code = load_template('kernel_iq.c')
735    user_code = [(f, open(f).read()) for f in model_sources(model_info)]
736
737    # Build initial sources
738    source = []
739    _add_source(source, *kernel_header)
740    for path, code in user_code:
741        _add_source(source, code, path)
742
743    if model_info.c_code:
744        _add_source(source, model_info.c_code, model_info.filename,
745                    lineno=model_info.lineno.get('c_code', 1))
746
747    # Make parameters for q, qx, qy so that we can use them in declarations
748    q, qx, qy, qab, qa, qb, qc \
749        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
750    # Generate form_volume function, etc. from body only
751    if isinstance(model_info.form_volume, str):
752        pars = partable.form_volume_parameters
753        source.append(_gen_fn(model_info, 'form_volume', pars))
754    if isinstance(model_info.Iq, str):
755        pars = [q] + partable.iq_parameters
756        source.append(_gen_fn(model_info, 'Iq', pars))
757    if isinstance(model_info.Iqxy, str):
758        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters
759        source.append(_gen_fn(model_info, 'Iqxy', pars))
760    if isinstance(model_info.Iqac, str):
761        pars = [qab, qc] + partable.iq_parameters
762        source.append(_gen_fn(model_info, 'Iqac', pars))
763    if isinstance(model_info.Iqabc, str):
764        pars = [qa, qb, qc] + partable.iq_parameters
765        source.append(_gen_fn(model_info, 'Iqabc', pars))
766
767    # What kind of 2D model do we need?  Is it consistent with the parameters?
768    xy_mode = find_xy_mode(source)
769    if xy_mode == 'qabc' and not partable.is_asymmetric:
770        raise ValueError("asymmetric oriented models need to define Iqabc")
771    elif xy_mode == 'qac' and partable.is_asymmetric:
772        raise ValueError("symmetric oriented models need to define Iqac")
773    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'):
774        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode)
775    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'):
776        if xy_mode == 'qxy':
777            logger.warn("oriented shapes should define Iqac or Iqabc")
778        else:
779            raise ValueError("Expected function Iqac or Iqabc for oriented shape")
780
781    # Define the parameter table
782    lineno = getframeinfo(currentframe()).lineno + 2
783    source.append('#line %d "sasmodels/generate.py"'%lineno)
784    #source.append('introduce breakage in generate to test lineno reporting')
785    source.append("#define PARAMETER_TABLE \\")
786    source.append("\\\n".join(p.as_definition()
787                              for p in partable.kernel_parameters))
788
789    # Define the function calls
790    if partable.form_volume_parameters:
791        refs = _call_pars("_v.", partable.form_volume_parameters)
792        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs))
793    else:
794        # Model doesn't have volume.  We could make the kernel run a little
795        # faster by not using/transferring the volume normalizations, but
796        # the ifdef's reduce readability more than is worthwhile.
797        call_volume = "#define CALL_VOLUME(v) 1.0"
798    source.append(call_volume)
799
800    model_refs = _call_pars("_v.", partable.iq_parameters)
801    pars = ",".join(["_q"] + model_refs)
802    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
803    if xy_mode == 'qabc':
804        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
805        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars
806        clear_iqxy = "#undef CALL_IQ_ABC"
807    elif xy_mode == 'qac':
808        pars = ",".join(["_qa", "_qc"] + model_refs)
809        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
810        clear_iqxy = "#undef CALL_IQ_AC"
811    elif xy_mode == 'qa':
812        pars = ",".join(["_qa"] + model_refs)
813        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
814        clear_iqxy = "#undef CALL_IQ_A"
815    elif xy_mode == 'qxy':
816        orientation_refs = _call_pars("_v.", partable.orientation_parameters)
817        pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs)
818        call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars
819        clear_iqxy = "#undef CALL_IQ_XY"
820        if partable.orientation_parameters:
821            call_iqxy += "\n#define HAVE_THETA"
822            clear_iqxy += "\n#undef HAVE_THETA"
823        if partable.is_asymmetric:
824            call_iqxy += "\n#define HAVE_PSI"
825            clear_iqxy += "\n#undef HAVE_PSI"
826
827
828    magpars = [k-2 for k, p in enumerate(partable.call_parameters)
829               if p.type == 'sld']
830
831    # Fill in definitions for numbers of parameters
832    source.append("#define MAX_PD %s"%partable.max_pd)
833    source.append("#define NUM_PARS %d"%partable.npars)
834    source.append("#define NUM_VALUES %d" % partable.nvalues)
835    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
836    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
837    source.append("#define PROJECTION %d"%PROJECTION)
838
839    # TODO: allow mixed python/opencl kernels?
840
841    ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
842    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
843    result = {
844        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
845        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
846    }
847
848    return result
849
850
851def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name):
852    # type: ([str,str], str, str, str) -> List[str]
853    code = kernel[0]
854    path = kernel[1].replace('\\', '\\\\')
855    iq = [
856        # define the Iq kernel
857        "#define KERNEL_NAME %s_Iq" % name,
858        call_iq,
859        '#line 1 "%s Iq"' % path,
860        code,
861        "#undef CALL_IQ",
862        "#undef KERNEL_NAME",
863        ]
864
865    iqxy = [
866        # define the Iqxy kernel from the same source with different #defines
867        "#define KERNEL_NAME %s_Iqxy" % name,
868        call_iqxy,
869        '#line 1 "%s Iqxy"' % path,
870        code,
871        clear_iqxy,
872        "#undef KERNEL_NAME",
873    ]
874
875    imagnetic = [
876        # define the Imagnetic kernel
877        "#define KERNEL_NAME %s_Imagnetic" % name,
878        "#define MAGNETIC 1",
879        call_iqxy,
880        '#line 1 "%s Imagnetic"' % path,
881        code,
882        clear_iqxy,
883        "#undef MAGNETIC",
884        "#undef KERNEL_NAME",
885    ]
886
887    return iq, iqxy, imagnetic
888
889
890def load_kernel_module(model_name):
891    # type: (str) -> module
892    """
893    Return the kernel module named in *model_name*.
894
895    If the name ends in *.py* then load it as a custom model using
896    :func:`sasmodels.custom.load_custom_kernel_module`, otherwise
897    load it from :mod:`sasmodels.models`.
898    """
899    if model_name.endswith('.py'):
900        kernel_module = load_custom_kernel_module(model_name)
901    else:
902        from sasmodels import models
903        __import__('sasmodels.models.'+model_name)
904        kernel_module = getattr(models, model_name, None)
905    return kernel_module
906
907
908section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z'
909                            % re.escape(string.punctuation))
910def _convert_section_titles_to_boldface(lines):
911    # type: (Sequence[str]) -> Iterator[str]
912    """
913    Do the actual work of identifying and converting section headings.
914    """
915    prior = None
916    for line in lines:
917        if prior is None:
918            prior = line
919        elif section_marker.match(line):
920            if len(line) >= len(prior):
921                yield "".join(("**", prior, "**"))
922                prior = None
923            else:
924                yield prior
925                prior = line
926        else:
927            yield prior
928            prior = line
929    if prior is not None:
930        yield prior
931
932
933def convert_section_titles_to_boldface(s):
934    # type: (str) -> str
935    """
936    Use explicit bold-face rather than section headings so that the table of
937    contents is not polluted with section names from the model documentation.
938
939    Sections are identified as the title line followed by a line of punctuation
940    at least as long as the title line.
941    """
942    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
943
944
945def make_doc(model_info):
946    # type: (ModelInfo) -> str
947    """
948    Return the documentation for the model.
949    """
950    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
951    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
952    docs = model_info.docs if model_info.docs is not None else ""
953    docs = convert_section_titles_to_boldface(docs)
954    pars = make_partable(model_info.parameters.COMMON
955                         + model_info.parameters.kernel_parameters)
956    subst = dict(id=model_info.id.replace('_', '-'),
957                 name=model_info.name,
958                 title=model_info.title,
959                 parameters=pars,
960                 returns=Sq_units if model_info.structure_factor else Iq_units,
961                 docs=docs)
962    return DOC_HEADER % subst
963
964
965# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog
966RST_PROLOG = r"""\
967.. |Ang| unicode:: U+212B
968.. |Ang^-1| replace:: |Ang|\ :sup:`-1`
969.. |Ang^2| replace:: |Ang|\ :sup:`2`
970.. |Ang^-2| replace:: |Ang|\ :sup:`-2`
971.. |1e-6Ang^-2| replace:: 10\ :sup:`-6`\ |Ang|\ :sup:`-2`
972.. |Ang^3| replace:: |Ang|\ :sup:`3`
973.. |Ang^-3| replace:: |Ang|\ :sup:`-3`
974.. |Ang^-4| replace:: |Ang|\ :sup:`-4`
975.. |cm^-1| replace:: cm\ :sup:`-1`
976.. |cm^2| replace:: cm\ :sup:`2`
977.. |cm^-2| replace:: cm\ :sup:`-2`
978.. |cm^3| replace:: cm\ :sup:`3`
979.. |1e15cm^3| replace:: 10\ :sup:`15`\ cm\ :sup:`3`
980.. |cm^-3| replace:: cm\ :sup:`-3`
981.. |sr^-1| replace:: sr\ :sup:`-1`
982
983.. |cdot| unicode:: U+00B7
984.. |deg| unicode:: U+00B0
985.. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3`
986.. |mg/m^2| replace:: mg\ |cdot|\ m\ :sup:`-2`
987.. |fm^2| replace:: fm\ :sup:`2`
988.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1`
989"""
990
991# TODO: make a better fake reference role
992RST_ROLES = """\
993.. role:: ref
994
995.. role:: numref
996
997"""
998
999def make_html(model_info):
1000    # type: (ModelInfo) -> str
1001    """
1002    Convert model docs directly to html.
1003    """
1004    from . import rst2html
1005
1006    rst = make_doc(model_info)
1007    return rst2html.rst2html("".join((RST_ROLES, RST_PROLOG, rst)))
1008
1009def view_html(model_name):
1010    # type: (str) -> None
1011    """
1012    Load the model definition and view its help.
1013    """
1014    from . import modelinfo
1015    kernel_module = load_kernel_module(model_name)
1016    info = modelinfo.make_model_info(kernel_module)
1017    view_html_from_info(info)
1018
1019def view_html_from_info(info):
1020    # type: (ModelInfo) -> None
1021    """
1022    View the help for a loaded model definition.
1023    """
1024    from . import rst2html
1025    url = "file://"+dirname(info.filename)+"/"
1026    rst2html.view_html(make_html(info), url=url)
1027
1028def demo_time():
1029    # type: () -> None
1030    """
1031    Show how long it takes to process a model.
1032    """
1033    import datetime
1034    from .modelinfo import make_model_info
1035    from .models import cylinder
1036
1037    tic = datetime.datetime.now()
1038    make_source(make_model_info(cylinder))
1039    toc = (datetime.datetime.now() - tic).total_seconds()
1040    print("time: %g"%toc)
1041
1042
1043def main():
1044    # type: () -> None
1045    """
1046    Program which prints the source produced by the model.
1047    """
1048    from .modelinfo import make_model_info
1049
1050    if len(sys.argv) <= 1:
1051        print("usage: python -m sasmodels.generate modelname")
1052    else:
1053        name = sys.argv[1]
1054        kernel_module = load_kernel_module(name)
1055        model_info = make_model_info(kernel_module)
1056        source = make_source(model_info)
1057        print(source['dll'])
1058
1059
1060if __name__ == "__main__":
1061    main()
Note: See TracBrowser for help on using the repository browser.