source: sasmodels/sasmodels/generate.py @ e2592f0

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

remove scale+background from S(q) parameter tables

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