source: sasmodels/sasmodels/generate.py @ 01c8d9e

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 01c8d9e was 01c8d9e, checked in by Suczewski <ges3@…>, 13 months ago

beta approximation, first pass

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