source: sasmodels/sasmodels/generate.py @ 9ee2756

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

simplify kernel wrapper code and combine OpenCL with DLL in one file

  • Property mode set to 100644
File size: 33.6 KB
Line 
1"""
2SAS model constructor.
3
4Small angle scattering models are defined by a set of kernel functions:
5
6    *Iq(q, p1, p2, ...)* returns the scattering at q for a form with
7    particular dimensions averaged over all orientations.
8
9    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form
10    with particular dimensions for a single orientation.
11
12    *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the
13    polarized neutron spin states (up-up, up-down, down-up, down-down) for
14    a form with particular dimensions for a single orientation.
15
16    *form_volume(p1, p2, ...)* returns the volume of the form with particular
17    dimension, or 1.0 if no volume normalization is required.
18
19    *ER(p1, p2, ...)* returns the effective radius of the form with
20    particular dimensions.
21
22    *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms.
23
24    #define INVALID(v) (expr)  returns False if v.parameter is invalid
25    for some parameter or other (e.g., v.bell_radius < v.radius).  If
26    necessary, the expression can call a function.
27
28These functions are defined in a kernel module .py script and an associated
29set of .c files.  The model constructor will use them to create models with
30polydispersity across volume and orientation parameters, and provide
31scale and background parameters for each model.
32
33*Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99
34functions written for OpenCL.  All functions need prototype declarations
35even if the are defined before they are used.  OpenCL does not support
36*#include* preprocessor directives, so instead the list of includes needs
37to be given as part of the metadata in the kernel module definition.
38The included files should be listed using a path relative to the kernel
39module, or if using "lib/file.c" if it is one of the standard includes
40provided with the sasmodels source.  The includes need to be listed in
41order so that functions are defined before they are used.
42
43Floating point values should be declared as *double*.  For single precision
44calculations, *double* will be replaced by *float*.  The single precision
45conversion will also tag floating point constants with "f" to make them
46single precision constants.  When using integral values in floating point
47expressions, they should be expressed as floating point values by including
48a decimal point.  This includes 0., 1. and 2.
49
50OpenCL has a *sincos* function which can improve performance when both
51the *sin* and *cos* values are needed for a particular argument.  Since
52this function does not exist in C99, all use of *sincos* should be
53replaced by the macro *SINCOS(value, sn, cn)* where *sn* and *cn* are
54previously declared *double* variables.  When compiled for systems without
55OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.   If *value* is
56an expression, it will appear twice in this case; whether or not it will be
57evaluated twice depends on the quality of the compiler.
58
59If the input parameters are invalid, the scattering calculator should
60return a negative number. Particularly with polydispersity, there are
61some sets of shape parameters which lead to nonsensical forms, such
62as a capped cylinder where the cap radius is smaller than the
63cylinder radius.  The polydispersity calculation will ignore these points,
64effectively chopping the parameter weight distributions at the boundary
65of the infeasible region.  The resulting scattering will be set to
66background.  This will work correctly even when polydispersity is off.
67
68*ER* and *VR* are python functions which operate on parameter vectors.
69The constructor code will generate the necessary vectors for computing
70them with the desired polydispersity.
71The kernel module must set variables defining the kernel meta data:
72
73    *id* is an implicit variable formed from the filename.  It will be
74    a valid python identifier, and will be used as the reference into
75    the html documentation, with '_' replaced by '-'.
76
77    *name* is the model name as displayed to the user.  If it is missing,
78    it will be constructed from the id.
79
80    *title* is a short description of the model, suitable for a tool tip,
81    or a one line model summary in a table of models.
82
83    *description* is an extended description of the model to be displayed
84    while the model parameters are being edited.
85
86    *parameters* is the list of parameters.  Parameters in the kernel
87    functions must appear in the same order as they appear in the
88    parameters list.  Two additional parameters, *scale* and *background*
89    are added to the beginning of the parameter list.  They will show up
90    in the documentation as model parameters, but they are never sent to
91    the kernel functions.  Note that *effect_radius* and *volfraction*
92    must occur first in structure factor calculations.
93
94    *category* is the default category for the model.  The category is
95    two level structure, with the form "group:section", indicating where
96    in the manual the model will be located.  Models are alphabetical
97    within their section.
98
99    *source* is the list of C-99 source files that must be joined to
100    create the OpenCL kernel functions.  The files defining the functions
101    need to be listed before the files which use the functions.
102
103    *ER* is a python function defining the effective radius.  If it is
104    not present, the effective radius is 0.
105
106    *VR* is a python function defining the volume ratio.  If it is not
107    present, the volume ratio is 1.
108
109    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the
110    C source code for the body of the volume, Iq, and Iqxy functions
111    respectively.  These can also be defined in the last source file.
112
113    *Iq* and *Iqxy* also be instead be python functions defining the
114    kernel.  If they are marked as *Iq.vectorized = True* then the
115    kernel is passed the entire *q* vector at once, otherwise it is
116    passed values one *q* at a time.  The performance improvement of
117    this step is significant.
118
119    *demo* is a dictionary of parameter=value defining a set of
120    parameters to use by default when *compare* is called.  Any
121    parameter not set in *demo* gets the initial value from the
122    parameter list.  *demo* is mostly needed to set the default
123    polydispersity values for tests.
124
125A :class:`modelinfo.ModelInfo` structure is constructed from the kernel meta
126data and returned to the caller.
127
128The doc string at the start of the kernel module will be used to
129construct the model documentation web pages.  Embedded figures should
130appear in the subdirectory "img" beside the model definition, and tagged
131with the kernel module name to avoid collision with other models.  Some
132file systems are case-sensitive, so only use lower case characters for
133file names and extensions.
134
135Code follows the C99 standard with the following extensions and conditions::
136
137    M_PI_180 = pi/180
138    M_4PI_3 = 4pi/3
139    square(x) = x*x
140    cube(x) = x*x*x
141    sas_sinx_x(x) = sin(x)/x, with sin(0)/0 -> 1
142    all double precision constants must include the decimal point
143    all double declarations may be converted to half, float, or long double
144    FLOAT_SIZE is the number of bytes in the converted variables
145
146:func:`load_kernel_module` loads the model definition file and
147:func:`modelinfo.make_model_info` parses it. :func:`make_source`
148converts C-based model definitions to C source code, including the
149polydispersity integral.  :func:`model_sources` returns the list of
150source files the model depends on, and :func:`timestamp` returns
151the latest time stamp amongst the source files (so you can check if
152the model needs to be rebuilt).
153
154The function :func:`make_doc` extracts the doc string and adds the
155parameter table to the top.  *make_figure* in *sasmodels/doc/genmodel*
156creates the default figure for the model.  [These two sets of code
157should mignrate into docs.py so docs can be updated in one place].
158"""
159from __future__ import print_function
160
161# TODO: determine which functions are useful outside of generate
162#__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
163
164import sys
165from os.path import abspath, dirname, join as joinpath, exists, isdir, getmtime
166import re
167import string
168from zlib import crc32
169from inspect import currentframe, getframeinfo
170
171import numpy as np  # type: ignore
172
173from .modelinfo import Parameter
174from .custom import load_custom_kernel_module
175
176try:
177    from typing import Tuple, Sequence, Iterator, Dict
178    from .modelinfo import ModelInfo
179except ImportError:
180    pass
181
182def get_data_path(external_dir, target_file):
183    path = abspath(dirname(__file__))
184    if exists(joinpath(path, target_file)):
185        return path
186
187    # check next to exe/zip file
188    exepath = dirname(sys.executable)
189    path = joinpath(exepath, external_dir)
190    if exists(joinpath(path, target_file)):
191        return path
192
193    # check in py2app Contents/Resources
194    path = joinpath(exepath, '..', 'Resources', external_dir)
195    if exists(joinpath(path, target_file)):
196        return abspath(path)
197
198    raise RuntimeError('Could not find '+joinpath(external_dir, target_file))
199
200EXTERNAL_DIR = 'sasmodels-data'
201DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_template.c')
202MODEL_PATH = joinpath(DATA_PATH, 'models')
203
204F16 = np.dtype('float16')
205F32 = np.dtype('float32')
206F64 = np.dtype('float64')
207try:  # CRUFT: older numpy does not support float128
208    F128 = np.dtype('float128')
209except TypeError:
210    F128 = None
211
212# Conversion from units defined in the parameter table for each model
213# to units displayed in the sphinx documentation.
214# This section associates the unit with the macro to use to produce the LaTex
215# code.  The macro itself needs to be defined in sasmodels/doc/rst_prolog.
216#
217# NOTE: there is an RST_PROLOG at the end of this file which is NOT
218# used for the bundled documentation. Still as long as we are defining the macros
219# in two places any new addition should define the macro in both places.
220RST_UNITS = {
221    "Ang": "|Ang|",
222    "1/Ang": "|Ang^-1|",
223    "1/Ang^2": "|Ang^-2|",
224    "Ang^3": "|Ang^3|",
225    "Ang^2": "|Ang^2|",
226    "1e15/cm^3": "|1e15cm^3|",
227    "Ang^3/mol": "|Ang^3|/mol",
228    "1e-6/Ang^2": "|1e-6Ang^-2|",
229    "degrees": "degree",
230    "1/cm": "|cm^-1|",
231    "Ang/cm": "|Ang*cm^-1|",
232    "g/cm^3": "|g/cm^3|",
233    "mg/m^2": "|mg/m^2|",
234    "": "None",
235    }
236
237# Headers for the parameters tables in th sphinx documentation
238PARTABLE_HEADERS = [
239    "Parameter",
240    "Description",
241    "Units",
242    "Default value",
243    ]
244
245# Minimum width for a default value (this is shorter than the column header
246# width, so will be ignored).
247PARTABLE_VALUE_WIDTH = 10
248
249# Documentation header for the module, giving the model name, its short
250# description and its parameter table.  The remainder of the doc comes
251# from the module docstring.
252DOC_HEADER = """.. _%(id)s:
253
254%(name)s
255=======================================================
256
257%(title)s
258
259%(parameters)s
260
261%(returns)s
262
263%(docs)s
264"""
265
266
267def format_units(units):
268    # type: (str) -> str
269    """
270    Convert units into ReStructured Text format.
271    """
272    return "string" if isinstance(units, list) else RST_UNITS.get(units, units)
273
274
275def make_partable(pars):
276    # type: (List[Parameter]) -> str
277    """
278    Generate the parameter table to include in the sphinx documentation.
279    """
280    column_widths = [
281        max(len(p.name) for p in pars),
282        max(len(p.description) for p in pars),
283        max(len(format_units(p.units)) for p in pars),
284        PARTABLE_VALUE_WIDTH,
285        ]
286    column_widths = [max(w, len(h))
287                     for w, h in zip(column_widths, PARTABLE_HEADERS)]
288
289    sep = " ".join("="*w for w in column_widths)
290    lines = [
291        sep,
292        " ".join("%-*s" % (w, h)
293                 for w, h in zip(column_widths, PARTABLE_HEADERS)),
294        sep,
295        ]
296    for p in pars:
297        lines.append(" ".join([
298            "%-*s" % (column_widths[0], p.name),
299            "%-*s" % (column_widths[1], p.description),
300            "%-*s" % (column_widths[2], format_units(p.units)),
301            "%*g" % (column_widths[3], p.default),
302            ]))
303    lines.append(sep)
304    return "\n".join(lines)
305
306
307def _search(search_path, filename):
308    # type: (List[str], str) -> str
309    """
310    Find *filename* in *search_path*.
311
312    Raises ValueError if file does not exist.
313    """
314    for path in search_path:
315        target = joinpath(path, filename)
316        if exists(target):
317            return target
318    raise ValueError("%r not found in %s" % (filename, search_path))
319
320
321def model_sources(model_info):
322    # type: (ModelInfo) -> List[str]
323    """
324    Return a list of the sources file paths for the module.
325    """
326    search_path = [dirname(model_info.filename), MODEL_PATH]
327    return [_search(search_path, f) for f in model_info.source]
328
329
330def dll_timestamp(model_info):
331    # type: (ModelInfo) -> int
332    """
333    Return a timestamp for the model corresponding to the most recently
334    changed file or dependency.
335    """
336    # TODO: fails DRY; templates appear two places.
337    model_templates = [joinpath(DATA_PATH, filename)
338                       for filename in ('kernel_header.c', 'kernel_iq.c')]
339    source_files = (model_sources(model_info)
340                    + model_templates
341                    + [model_info.filename])
342    # Note: file may not exist when it is a standard model from library.zip
343    times = [getmtime(f) for f in source_files if exists(f)]
344    newest = max(times) if times else 0
345    return newest
346
347def ocl_timestamp(model_info):
348    # type: (ModelInfo) -> int
349    """
350    Return a timestamp for the model corresponding to the most recently
351    changed file or dependency.
352
353    Note that this does not look at the time stamps for the OpenCL header
354    information since that need not trigger a recompile of the DLL.
355    """
356    # TODO: fails DRY; templates appear two places.
357    model_templates = [joinpath(DATA_PATH, filename)
358                       for filename in ('kernel_header.c', 'kernel_iq.cl')]
359    source_files = (model_sources(model_info)
360                    + model_templates
361                    + [model_info.filename])
362    # Note: file may not exist when it is a standard model from library.zip
363    times = [getmtime(f) for f in source_files if exists(f)]
364    newest = max(times) if times else 0
365    return newest
366
367def tag_source(source):
368    # type: (str) -> str
369    """
370    Return a unique tag for the source code.
371    """
372    # Note: need 0xffffffff&val to force an unsigned 32-bit number
373    return "%08X"%(0xffffffff&crc32(source))
374
375def convert_type(source, dtype):
376    # type: (str, np.dtype) -> str
377    """
378    Convert code from double precision to the desired type.
379
380    Floating point constants are tagged with 'f' for single precision or 'L'
381    for long double precision.
382    """
383    source = _fix_tgmath_int(source)
384    if dtype == F16:
385        fbytes = 2
386        source = _convert_type(source, "half", "f")
387    elif dtype == F32:
388        fbytes = 4
389        source = _convert_type(source, "float", "f")
390    elif dtype == F64:
391        fbytes = 8
392        # no need to convert if it is already double
393    elif dtype == F128:
394        fbytes = 16
395        source = _convert_type(source, "long double", "L")
396    else:
397        raise ValueError("Unexpected dtype in source conversion: %s" % dtype)
398    return ("#define FLOAT_SIZE %d\n" % fbytes)+source
399
400
401def _convert_type(source, type_name, constant_flag):
402    # type: (str, str, str) -> str
403    """
404    Replace 'double' with *type_name* in *source*, tagging floating point
405    constants with *constant_flag*.
406    """
407    # Convert double keyword to float/long double/half.
408    # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16.
409    # Assume complex numbers are represented as cdouble which is typedef'd
410    # to double2.
411    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
412                    r'\1%s\2'%type_name, source)
413    source = _tag_float(source, constant_flag)
414    return source
415
416TGMATH_INT_RE = re.compile(r"""
417(?: # Non-capturing match; not lookbehind since pattern length is variable
418  \b              # word boundary
419   # various math functions
420  (a?(sin|cos|tan)h? | atan2
421   | erfc? | tgamma
422   | exp(2|10|m1)? | log(2|10|1p)? | pow[nr]? | sqrt | rsqrt | rootn
423   | fabs | fmax | fmin
424   )
425  \s*[(]\s*       # open parenthesis
426)
427[+-]?(0|[1-9]\d*) # integer
428(?=               # lookahead match: don't want to move from end of int
429  \s*[,)]         # comma or close parenthesis for end of argument
430)                 # end lookahead
431""", re.VERBOSE)
432def _fix_tgmath_int(source):
433    # type: (str) -> str
434    """
435    Replace f(integer) with f(integer.) for sin, cos, pow, etc.
436
437    OS X OpenCL complains that it can't resolve the type generic calls to
438    the standard math functions when they are called with integer constants,
439    but this does not happen with the Windows Intel driver for example.
440    To avoid confusion on the matrix marketplace, automatically promote
441    integers to floats if we recognize them in the source.
442
443    The specific functions we look for are:
444
445        trigonometric: sin, asin, sinh, asinh, etc., and atan2
446        exponential:   exp, exp2, exp10, expm1, log, log2, log10, logp1
447        power:         pow, pown, powr, sqrt, rsqrt, rootn
448        special:       erf, erfc, tgamma
449        float:         fabs, fmin, fmax
450
451    Note that we don't convert the second argument of dual argument
452    functions: atan2, fmax, fmin, pow, powr.  This could potentially
453    be a problem for pow(x, 2), but that case seems to work without change.
454    """
455    out = TGMATH_INT_RE.sub(r'\g<0>.', source)
456    return out
457
458
459# Floating point regular expression
460#
461# Define parts:
462#
463#    E = [eE][+-]?\d+    : Exponent
464#    P = [.]             : Decimal separator
465#    N = [1-9]\d*        : Natural number, no leading zeros
466#    Z = 0               : Zero
467#    F = \d+             : Fractional number, maybe leading zeros
468#    F? = \d*            : Optional fractional number
469#
470# We want to reject bare natural numbers and bare decimal points, so we
471# need to tediously outline the cases where we have either a fraction or
472# an exponent:
473#
474#   ( ZP | ZPF | ZE | ZPE | ZPFE | NP | NPF | NE | NPE | NPFE | PF | PFE )
475#
476#
477# We can then join cases by making parts optional.  The following are
478# some ways to do this:
479#
480#   ( (Z|N)(P|PF|E|PE|PFE) | PFE? )                   # Split on lead
481#     => ( (Z|N)(PF?|(PF?)?E) | PFE? )
482#   ( ((Z|N)PF?|PF)E? | (Z|N)E)                       # Split on point
483#   ( (ZP|ZPF|NP|NPF|PF) | (Z|ZP|ZPF|N|NP|NPF|PF)E )  # Split on E
484#     => ( ((Z|N)PF?|PF) | ((Z|N)(PF?)? | PF) E )
485FLOAT_RE = re.compile(r"""
486    (?<!\w)  # use negative lookbehind since '.' confuses \b test
487    # use split on lead to match float ( (Z|N)(PF?|(PF?)?E) | PFE? )
488    ( ( 0 | [1-9]\d* )                     # ( ( Z | N )
489      ([.]\d* | ([.]\d*)? [eE][+-]?\d+ )   #   (PF? | (PF?)? E )
490    | [.]\d+ ([eE][+-]?\d+)?               # | PF (E)?
491    )                                      # )
492    (?!\w)  # use negative lookahead since '.' confuses \b test
493    """, re.VERBOSE)
494def _tag_float(source, constant_flag):
495    # Convert floating point constants to single by adding 'f' to the end,
496    # or long double with an 'L' suffix.  OS/X complains if you don't do this.
497    out = FLOAT_RE.sub(r'\g<0>%s'%constant_flag, source)
498    #print("in",repr(source),"out",repr(out), constant_flag)
499    return out
500
501def test_tag_float():
502    """check that floating point constants are properly identified and tagged with 'f'"""
503
504    cases = """
505ZP  : 0.
506ZPF : 0.0,0.01,0.1
507Z  E: 0e+001
508ZP E: 0.E0
509ZPFE: 0.13e-031
510NP  : 1., 12.
511NPF : 1.0001, 1.1, 1.0
512N  E: 1e0, 37E-080
513NP E: 1.e0, 37.E-080
514NPFE: 845.017e+22
515 PF : .1, .0, .0100
516 PFE: .6e+9, .82E-004
517# isolated cases
5180.
5191e0
5200.13e-013
521# untouched
522struct3.e3, 03.05.67, 37
523# expressions
5243.75+-1.6e-7-27+13.2
525a3.e2 - 0.
5264*atan(1)
5274.*atan(1.)
528"""
529
530    output = """
531ZP  : 0.f
532ZPF : 0.0f,0.01f,0.1f
533Z  E: 0e+001f
534ZP E: 0.E0f
535ZPFE: 0.13e-031f
536NP  : 1.f, 12.f
537NPF : 1.0001f, 1.1f, 1.0f
538N  E: 1e0f, 37E-080f
539NP E: 1.e0f, 37.E-080f
540NPFE: 845.017e+22f
541 PF : .1f, .0f, .0100f
542 PFE: .6e+9f, .82E-004f
543# isolated cases
5440.f
5451e0f
5460.13e-013f
547# untouched
548struct3.e3, 03.05.67, 37
549# expressions
5503.75f+-1.6e-7f-27+13.2f
551a3.e2 - 0.f
5524*atan(1)
5534.f*atan(1.f)
554"""
555
556    for case_in, case_out in zip(cases.split('\n'), output.split('\n')):
557        out = _tag_float(case_in, 'f')
558        assert case_out == out, "%r => %r"%(case_in, out)
559
560
561def kernel_name(model_info, variant):
562    # type: (ModelInfo, str) -> str
563    """
564    Name of the exported kernel symbol.
565
566    *variant* is "Iq", "Iqxy" or "Imagnetic".
567    """
568    return model_info.name + "_" + variant
569
570
571def indent(s, depth):
572    # type: (str, int) -> str
573    """
574    Indent a string of text with *depth* additional spaces on each line.
575    """
576    spaces = " "*depth
577    sep = "\n" + spaces
578    return spaces + sep.join(s.split("\n"))
579
580
581_template_cache = {}  # type: Dict[str, Tuple[int, str, str]]
582def load_template(filename):
583    # type: (str) -> str
584    path = joinpath(DATA_PATH, filename)
585    mtime = getmtime(path)
586    if filename not in _template_cache or mtime > _template_cache[filename][0]:
587        with open(path) as fid:
588            _template_cache[filename] = (mtime, fid.read(), path)
589    return _template_cache[filename][1], path
590
591
592_FN_TEMPLATE = """\
593double %(name)s(%(pars)s);
594double %(name)s(%(pars)s) {
595#line %(line)d "%(filename)s"
596    %(body)s
597}
598
599"""
600def _gen_fn(name, pars, body, filename, line):
601    # type: (str, List[Parameter], str, str, int) -> str
602    """
603    Generate a function given pars and body.
604
605    Returns the following string::
606
607         double fn(double a, double b, ...);
608         double fn(double a, double b, ...) {
609             ....
610         }
611    """
612    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
613    return _FN_TEMPLATE % {
614        'name': name, 'pars': par_decl, 'body': body,
615        'filename': filename.replace('\\', '\\\\'), 'line': line,
616    }
617
618
619def _call_pars(prefix, pars):
620    # type: (str, List[Parameter]) -> List[str]
621    """
622    Return a list of *prefix+parameter* from parameter items.
623
624    *prefix* should be "v." if v is a struct.
625    """
626    return [p.as_call_reference(prefix) for p in pars]
627
628
629# type in IQXY pattern could be single, float, double, long double, ...
630_IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)",
631                           flags=re.MULTILINE)
632def _have_Iqxy(sources):
633    # type: (List[str]) -> bool
634    """
635    Return true if any file defines Iqxy.
636
637    Note this is not a C parser, and so can be easily confused by
638    non-standard syntax.  Also, it will incorrectly identify the following
639    as having Iqxy::
640
641        /*
642        double Iqxy(qx, qy, ...) { ... fill this in later ... }
643        */
644
645    If you want to comment out an Iqxy function, use // on the front of the
646    line instead.
647    """
648    for path, code in sources:
649        if _IQXY_PATTERN.search(code):
650            return True
651    else:
652        return False
653
654
655def _add_source(source, code, path):
656    """
657    Add a file to the list of source code chunks, tagged with path and line.
658    """
659    path = path.replace('\\', '\\\\')
660    source.append('#line 1 "%s"' % path)
661    source.append(code)
662
663def make_source(model_info):
664    # type: (ModelInfo) -> Dict[str, str]
665    """
666    Generate the OpenCL/ctypes kernel from the module info.
667
668    Uses source files found in the given search path.  Returns None if this
669    is a pure python model, with no C source components.
670    """
671    if callable(model_info.Iq):
672        raise ValueError("can't compile python model")
673        #return None
674
675    # TODO: need something other than volume to indicate dispersion parameters
676    # No volume normalization despite having a volume parameter.
677    # Thickness is labelled a volume in order to trigger polydispersity.
678    # May want a separate dispersion flag, or perhaps a separate category for
679    # disperse, but not volume.  Volume parameters also use relative values
680    # for the distribution rather than the absolute values used by angular
681    # dispersion.  Need to be careful that necessary parameters are available
682    # for computing volume even if we allow non-disperse volume parameters.
683
684    partable = model_info.parameters
685
686    # Load templates and user code
687    kernel_header = load_template('kernel_header.c')
688    kernel_code = load_template('kernel_iq.c')
689    user_code = [(f, open(f).read()) for f in model_sources(model_info)]
690
691    # What kind of 2D model do we need?
692    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str)
693               else 'qac' if not partable.is_asymmetric
694               else 'qabc')
695
696    # Build initial sources
697    source = []
698    _add_source(source, *kernel_header)
699    for path, code in user_code:
700        _add_source(source, code, path)
701
702    # Make parameters for q, qx, qy so that we can use them in declarations
703    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')]
704    # Generate form_volume function, etc. from body only
705    if isinstance(model_info.form_volume, str):
706        pars = partable.form_volume_parameters
707        source.append(_gen_fn('form_volume', pars, model_info.form_volume,
708                              model_info.filename, model_info._form_volume_line))
709    if isinstance(model_info.Iq, str):
710        pars = [q] + partable.iq_parameters
711        source.append(_gen_fn('Iq', pars, model_info.Iq,
712                              model_info.filename, model_info._Iq_line))
713    if isinstance(model_info.Iqxy, str):
714        pars = [qx, qy] + partable.iqxy_parameters
715        source.append(_gen_fn('Iqxy', pars, model_info.Iqxy,
716                              model_info.filename, model_info._Iqxy_line))
717
718    # Define the parameter table
719    lineno = getframeinfo(currentframe()).lineno + 2
720    source.append('#line %d "sasmodels/generate.py"'%lineno)
721    #source.append('introduce breakage in generate to test lineno reporting')
722    source.append("#define PARAMETER_TABLE \\")
723    source.append("\\\n".join(p.as_definition()
724                              for p in partable.kernel_parameters))
725
726    # Define the function calls
727    if partable.form_volume_parameters:
728        refs = _call_pars("_v.", partable.form_volume_parameters)
729        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs))
730    else:
731        # Model doesn't have volume.  We could make the kernel run a little
732        # faster by not using/transferring the volume normalizations, but
733        # the ifdef's reduce readability more than is worthwhile.
734        call_volume = "#define CALL_VOLUME(v) 1.0"
735    source.append(call_volume)
736
737    model_refs = _call_pars("_v.", partable.iq_parameters)
738    pars = ",".join(["_q"] + model_refs)
739    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
740    if xy_mode == 'qabc':
741        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
742        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars
743        clear_iqxy = "#undef CALL_IQ_ABC"
744    elif xy_mode == 'qac':
745        pars = ",".join(["_qa", "_qc"] + model_refs)
746        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars
747        clear_iqxy = "#undef CALL_IQ_AC"
748    else:  # xy_mode == 'qa'
749        pars = ",".join(["_qa"] + model_refs)
750        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
751        clear_iqxy = "#undef CALL_IQ_A"
752
753    magpars = [k-2 for k, p in enumerate(partable.call_parameters)
754               if p.type == 'sld']
755
756    # Fill in definitions for numbers of parameters
757    source.append("#define MAX_PD %s"%partable.max_pd)
758    source.append("#define NUM_PARS %d"%partable.npars)
759    source.append("#define NUM_VALUES %d" % partable.nvalues)
760    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
761    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
762    for k, v in enumerate(magpars[:3]):
763        source.append("#define MAGNETIC_PAR%d %d"%(k+1, v))
764
765    # TODO: allow mixed python/opencl kernels?
766
767    ocl = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
768    dll = kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
769    result = {
770        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
771        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
772    }
773
774    return result
775
776
777def kernels(kernel, call_iq, call_iqxy, clear_iqxy, name):
778    # type: ([str,str], str, str, str) -> List[str]
779    code = kernel[0]
780    path = kernel[1].replace('\\', '\\\\')
781    iq = [
782        # define the Iq kernel
783        "#define KERNEL_NAME %s_Iq" % name,
784        call_iq,
785        '#line 1 "%s Iq"' % path,
786        code,
787        "#undef CALL_IQ",
788        "#undef KERNEL_NAME",
789        ]
790
791    iqxy = [
792        # define the Iqxy kernel from the same source with different #defines
793        "#define KERNEL_NAME %s_Iqxy" % name,
794        call_iqxy,
795        '#line 1 "%s Iqxy"' % path,
796        code,
797        clear_iqxy,
798        "#undef KERNEL_NAME",
799    ]
800
801    imagnetic = [
802        # define the Imagnetic kernel
803        "#define KERNEL_NAME %s_Imagnetic" % name,
804        "#define MAGNETIC 1",
805        call_iqxy,
806        '#line 1 "%s Imagnetic"' % path,
807        code,
808        clear_iqxy,
809        "#undef MAGNETIC",
810        "#undef KERNEL_NAME",
811    ]
812
813    return iq, iqxy, imagnetic
814
815
816def load_kernel_module(model_name):
817    # type: (str) -> module
818    """
819    Return the kernel module named in *model_name*.
820
821    If the name ends in *.py* then load it as a custom model using
822    :func:`sasmodels.custom.load_custom_kernel_module`, otherwise
823    load it from :mod:`sasmodels.models`.
824    """
825    if model_name.endswith('.py'):
826        kernel_module = load_custom_kernel_module(model_name)
827    else:
828        from sasmodels import models
829        __import__('sasmodels.models.'+model_name)
830        kernel_module = getattr(models, model_name, None)
831    return kernel_module
832
833
834section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z'
835                            % re.escape(string.punctuation))
836def _convert_section_titles_to_boldface(lines):
837    # type: (Sequence[str]) -> Iterator[str]
838    """
839    Do the actual work of identifying and converting section headings.
840    """
841    prior = None
842    for line in lines:
843        if prior is None:
844            prior = line
845        elif section_marker.match(line):
846            if len(line) >= len(prior):
847                yield "".join(("**", prior, "**"))
848                prior = None
849            else:
850                yield prior
851                prior = line
852        else:
853            yield prior
854            prior = line
855    if prior is not None:
856        yield prior
857
858
859def convert_section_titles_to_boldface(s):
860    # type: (str) -> str
861    """
862    Use explicit bold-face rather than section headings so that the table of
863    contents is not polluted with section names from the model documentation.
864
865    Sections are identified as the title line followed by a line of punctuation
866    at least as long as the title line.
867    """
868    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
869
870
871def make_doc(model_info):
872    # type: (ModelInfo) -> str
873    """
874    Return the documentation for the model.
875    """
876    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
877    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
878    docs = model_info.docs if model_info.docs is not None else ""
879    docs = convert_section_titles_to_boldface(docs)
880    pars = make_partable(model_info.parameters.COMMON
881                         + model_info.parameters.kernel_parameters)
882    subst = dict(id=model_info.id.replace('_', '-'),
883                 name=model_info.name,
884                 title=model_info.title,
885                 parameters=pars,
886                 returns=Sq_units if model_info.structure_factor else Iq_units,
887                 docs=docs)
888    return DOC_HEADER % subst
889
890
891# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog
892RST_PROLOG = r"""\
893.. |Ang| unicode:: U+212B
894.. |Ang^-1| replace:: |Ang|\ :sup:`-1`
895.. |Ang^2| replace:: |Ang|\ :sup:`2`
896.. |Ang^-2| replace:: |Ang|\ :sup:`-2`
897.. |1e-6Ang^-2| replace:: 10\ :sup:`-6`\ |Ang|\ :sup:`-2`
898.. |Ang^3| replace:: |Ang|\ :sup:`3`
899.. |Ang^-3| replace:: |Ang|\ :sup:`-3`
900.. |Ang^-4| replace:: |Ang|\ :sup:`-4`
901.. |cm^-1| replace:: cm\ :sup:`-1`
902.. |cm^2| replace:: cm\ :sup:`2`
903.. |cm^-2| replace:: cm\ :sup:`-2`
904.. |cm^3| replace:: cm\ :sup:`3`
905.. |1e15cm^3| replace:: 10\ :sup:`15`\ cm\ :sup:`3`
906.. |cm^-3| replace:: cm\ :sup:`-3`
907.. |sr^-1| replace:: sr\ :sup:`-1`
908
909.. |cdot| unicode:: U+00B7
910.. |deg| unicode:: U+00B0
911.. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3`
912.. |mg/m^2| replace:: mg\ |cdot|\ m\ :sup:`-2`
913.. |fm^2| replace:: fm\ :sup:`2`
914.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1`
915"""
916
917# TODO: make a better fake reference role
918RST_ROLES = """\
919.. role:: ref
920
921.. role:: numref
922
923"""
924
925def make_html(model_info):
926    """
927    Convert model docs directly to html.
928    """
929    from . import rst2html
930
931    rst = make_doc(model_info)
932    return rst2html.rst2html("".join((RST_ROLES, RST_PROLOG, rst)))
933
934def view_html(model_name):
935    from . import modelinfo
936    kernel_module = load_kernel_module(model_name)
937    info = modelinfo.make_model_info(kernel_module)
938    view_html_from_info(info)
939
940def view_html_from_info(info):
941    from . import rst2html
942    url = "file://"+dirname(info.filename)+"/"
943    rst2html.view_html(make_html(info), url=url)
944
945def demo_time():
946    # type: () -> None
947    """
948    Show how long it takes to process a model.
949    """
950    import datetime
951    from .modelinfo import make_model_info
952    from .models import cylinder
953
954    tic = datetime.datetime.now()
955    make_source(make_model_info(cylinder))
956    toc = (datetime.datetime.now() - tic).total_seconds()
957    print("time: %g"%toc)
958
959
960def main():
961    # type: () -> None
962    """
963    Program which prints the source produced by the model.
964    """
965    import sys
966    from .modelinfo import make_model_info
967
968    if len(sys.argv) <= 1:
969        print("usage: python -m sasmodels.generate modelname")
970    else:
971        name = sys.argv[1]
972        kernel_module = load_kernel_module(name)
973        model_info = make_model_info(kernel_module)
974        source = make_source(model_info)
975        print(source['dll'])
976
977
978if __name__ == "__main__":
979    main()
Note: See TracBrowser for help on using the repository browser.