source: sasmodels/sasmodels/generate.py @ 2d81cfe

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

lint

  • Property mode set to 100644
File size: 34.1 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, 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
176# pylint: disable=unused-import
177try:
178    from typing import Tuple, Sequence, Iterator, Dict
179    from .modelinfo import ModelInfo
180except ImportError:
181    pass
182# pylint: enable=unused-import
183
184# jitter projection to use in the kernel code.  See explore/jitter.py
185# for details.  To change it from a program, set generate.PROJECTION.
186PROJECTION = 1
187
188def get_data_path(external_dir, target_file):
189    path = abspath(dirname(__file__))
190    if exists(joinpath(path, target_file)):
191        return path
192
193    # check next to exe/zip file
194    exepath = dirname(sys.executable)
195    path = joinpath(exepath, external_dir)
196    if exists(joinpath(path, target_file)):
197        return path
198
199    # check in py2app Contents/Resources
200    path = joinpath(exepath, '..', 'Resources', external_dir)
201    if exists(joinpath(path, target_file)):
202        return abspath(path)
203
204    raise RuntimeError('Could not find '+joinpath(external_dir, target_file))
205
206EXTERNAL_DIR = 'sasmodels-data'
207DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_template.c')
208MODEL_PATH = joinpath(DATA_PATH, 'models')
209
210F16 = np.dtype('float16')
211F32 = np.dtype('float32')
212F64 = np.dtype('float64')
213try:  # CRUFT: older numpy does not support float128
214    F128 = np.dtype('float128')
215except TypeError:
216    F128 = None
217
218# Conversion from units defined in the parameter table for each model
219# to units displayed in the sphinx documentation.
220# This section associates the unit with the macro to use to produce the LaTex
221# code.  The macro itself needs to be defined in sasmodels/doc/rst_prolog.
222#
223# NOTE: there is an RST_PROLOG at the end of this file which is NOT
224# used for the bundled documentation. Still as long as we are defining the macros
225# in two places any new addition should define the macro in both places.
226RST_UNITS = {
227    "Ang": "|Ang|",
228    "1/Ang": "|Ang^-1|",
229    "1/Ang^2": "|Ang^-2|",
230    "Ang^3": "|Ang^3|",
231    "Ang^2": "|Ang^2|",
232    "1e15/cm^3": "|1e15cm^3|",
233    "Ang^3/mol": "|Ang^3|/mol",
234    "1e-6/Ang^2": "|1e-6Ang^-2|",
235    "degrees": "degree",
236    "1/cm": "|cm^-1|",
237    "Ang/cm": "|Ang*cm^-1|",
238    "g/cm^3": "|g/cm^3|",
239    "mg/m^2": "|mg/m^2|",
240    "": "None",
241    }
242
243# Headers for the parameters tables in th sphinx documentation
244PARTABLE_HEADERS = [
245    "Parameter",
246    "Description",
247    "Units",
248    "Default value",
249    ]
250
251# Minimum width for a default value (this is shorter than the column header
252# width, so will be ignored).
253PARTABLE_VALUE_WIDTH = 10
254
255# Documentation header for the module, giving the model name, its short
256# description and its parameter table.  The remainder of the doc comes
257# from the module docstring.
258DOC_HEADER = """.. _%(id)s:
259
260%(name)s
261=======================================================
262
263%(title)s
264
265%(parameters)s
266
267%(returns)s
268
269%(docs)s
270"""
271
272
273def format_units(units):
274    # type: (str) -> str
275    """
276    Convert units into ReStructured Text format.
277    """
278    return "string" if isinstance(units, list) else RST_UNITS.get(units, units)
279
280
281def make_partable(pars):
282    # type: (List[Parameter]) -> str
283    """
284    Generate the parameter table to include in the sphinx documentation.
285    """
286    column_widths = [
287        max(len(p.name) for p in pars),
288        max(len(p.description) for p in pars),
289        max(len(format_units(p.units)) for p in pars),
290        PARTABLE_VALUE_WIDTH,
291        ]
292    column_widths = [max(w, len(h))
293                     for w, h in zip(column_widths, PARTABLE_HEADERS)]
294
295    sep = " ".join("="*w for w in column_widths)
296    lines = [
297        sep,
298        " ".join("%-*s" % (w, h)
299                 for w, h in zip(column_widths, PARTABLE_HEADERS)),
300        sep,
301        ]
302    for p in pars:
303        lines.append(" ".join([
304            "%-*s" % (column_widths[0], p.name),
305            "%-*s" % (column_widths[1], p.description),
306            "%-*s" % (column_widths[2], format_units(p.units)),
307            "%*g" % (column_widths[3], p.default),
308            ]))
309    lines.append(sep)
310    return "\n".join(lines)
311
312
313def _search(search_path, filename):
314    # type: (List[str], str) -> str
315    """
316    Find *filename* in *search_path*.
317
318    Raises ValueError if file does not exist.
319    """
320    for path in search_path:
321        target = joinpath(path, filename)
322        if exists(target):
323            return target
324    raise ValueError("%r not found in %s" % (filename, search_path))
325
326
327def model_sources(model_info):
328    # type: (ModelInfo) -> List[str]
329    """
330    Return a list of the sources file paths for the module.
331    """
332    search_path = [dirname(model_info.filename), MODEL_PATH]
333    return [_search(search_path, f) for f in model_info.source]
334
335
336def dll_timestamp(model_info):
337    # type: (ModelInfo) -> int
338    """
339    Return a timestamp for the model corresponding to the most recently
340    changed file or dependency.
341    """
342    # TODO: fails DRY; templates appear two places.
343    model_templates = [joinpath(DATA_PATH, filename)
344                       for filename in ('kernel_header.c', 'kernel_iq.c')]
345    source_files = (model_sources(model_info)
346                    + model_templates
347                    + [model_info.filename])
348    # Note: file may not exist when it is a standard model from library.zip
349    times = [getmtime(f) for f in source_files if exists(f)]
350    newest = max(times) if times else 0
351    return newest
352
353def ocl_timestamp(model_info):
354    # type: (ModelInfo) -> int
355    """
356    Return a timestamp for the model corresponding to the most recently
357    changed file or dependency.
358
359    Note that this does not look at the time stamps for the OpenCL header
360    information since that need not trigger a recompile of the DLL.
361    """
362    # TODO: fails DRY; templates appear two places.
363    model_templates = [joinpath(DATA_PATH, filename)
364                       for filename in ('kernel_header.c', 'kernel_iq.cl')]
365    source_files = (model_sources(model_info)
366                    + model_templates
367                    + [model_info.filename])
368    # Note: file may not exist when it is a standard model from library.zip
369    times = [getmtime(f) for f in source_files if exists(f)]
370    newest = max(times) if times else 0
371    return newest
372
373def tag_source(source):
374    # type: (str) -> str
375    """
376    Return a unique tag for the source code.
377    """
378    # Note: need 0xffffffff&val to force an unsigned 32-bit number
379    try:
380        source = source.encode('utf8')
381    except AttributeError: # bytes has no encode attribute in python 3
382        pass
383    return "%08X"%(0xffffffff&crc32(source))
384
385def convert_type(source, dtype):
386    # type: (str, np.dtype) -> str
387    """
388    Convert code from double precision to the desired type.
389
390    Floating point constants are tagged with 'f' for single precision or 'L'
391    for long double precision.
392    """
393    source = _fix_tgmath_int(source)
394    if dtype == F16:
395        fbytes = 2
396        source = _convert_type(source, "half", "f")
397    elif dtype == F32:
398        fbytes = 4
399        source = _convert_type(source, "float", "f")
400    elif dtype == F64:
401        fbytes = 8
402        # no need to convert if it is already double
403    elif dtype == F128:
404        fbytes = 16
405        source = _convert_type(source, "long double", "L")
406    else:
407        raise ValueError("Unexpected dtype in source conversion: %s" % dtype)
408    return ("#define FLOAT_SIZE %d\n" % fbytes)+source
409
410
411def _convert_type(source, type_name, constant_flag):
412    # type: (str, str, str) -> str
413    """
414    Replace 'double' with *type_name* in *source*, tagging floating point
415    constants with *constant_flag*.
416    """
417    # Convert double keyword to float/long double/half.
418    # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16.
419    # Assume complex numbers are represented as cdouble which is typedef'd
420    # to double2.
421    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
422                    r'\1%s\2'%type_name, source)
423    source = _tag_float(source, constant_flag)
424    return source
425
426TGMATH_INT_RE = re.compile(r"""
427(?: # Non-capturing match; not lookbehind since pattern length is variable
428  \b              # word boundary
429   # various math functions
430  (a?(sin|cos|tan)h? | atan2
431   | erfc? | tgamma
432   | exp(2|10|m1)? | log(2|10|1p)? | pow[nr]? | sqrt | rsqrt | rootn
433   | fabs | fmax | fmin
434   )
435  \s*[(]\s*       # open parenthesis
436)
437[+-]?(0|[1-9]\d*) # integer
438(?=               # lookahead match: don't want to move from end of int
439  \s*[,)]         # comma or close parenthesis for end of argument
440)                 # end lookahead
441""", re.VERBOSE)
442def _fix_tgmath_int(source):
443    # type: (str) -> str
444    """
445    Replace f(integer) with f(integer.) for sin, cos, pow, etc.
446
447    OS X OpenCL complains that it can't resolve the type generic calls to
448    the standard math functions when they are called with integer constants,
449    but this does not happen with the Windows Intel driver for example.
450    To avoid confusion on the matrix marketplace, automatically promote
451    integers to floats if we recognize them in the source.
452
453    The specific functions we look for are:
454
455        trigonometric: sin, asin, sinh, asinh, etc., and atan2
456        exponential:   exp, exp2, exp10, expm1, log, log2, log10, logp1
457        power:         pow, pown, powr, sqrt, rsqrt, rootn
458        special:       erf, erfc, tgamma
459        float:         fabs, fmin, fmax
460
461    Note that we don't convert the second argument of dual argument
462    functions: atan2, fmax, fmin, pow, powr.  This could potentially
463    be a problem for pow(x, 2), but that case seems to work without change.
464    """
465    out = TGMATH_INT_RE.sub(r'\g<0>.', source)
466    return out
467
468
469# Floating point regular expression
470#
471# Define parts:
472#
473#    E = [eE][+-]?\d+    : Exponent
474#    P = [.]             : Decimal separator
475#    N = [1-9]\d*        : Natural number, no leading zeros
476#    Z = 0               : Zero
477#    F = \d+             : Fractional number, maybe leading zeros
478#    F? = \d*            : Optional fractional number
479#
480# We want to reject bare natural numbers and bare decimal points, so we
481# need to tediously outline the cases where we have either a fraction or
482# an exponent:
483#
484#   ( ZP | ZPF | ZE | ZPE | ZPFE | NP | NPF | NE | NPE | NPFE | PF | PFE )
485#
486#
487# We can then join cases by making parts optional.  The following are
488# some ways to do this:
489#
490#   ( (Z|N)(P|PF|E|PE|PFE) | PFE? )                   # Split on lead
491#     => ( (Z|N)(PF?|(PF?)?E) | PFE? )
492#   ( ((Z|N)PF?|PF)E? | (Z|N)E)                       # Split on point
493#   ( (ZP|ZPF|NP|NPF|PF) | (Z|ZP|ZPF|N|NP|NPF|PF)E )  # Split on E
494#     => ( ((Z|N)PF?|PF) | ((Z|N)(PF?)? | PF) E )
495FLOAT_RE = re.compile(r"""
496    (?<!\w)  # use negative lookbehind since '.' confuses \b test
497    # use split on lead to match float ( (Z|N)(PF?|(PF?)?E) | PFE? )
498    ( ( 0 | [1-9]\d* )                     # ( ( Z | N )
499      ([.]\d* | ([.]\d*)? [eE][+-]?\d+ )   #   (PF? | (PF?)? E )
500    | [.]\d+ ([eE][+-]?\d+)?               # | PF (E)?
501    )                                      # )
502    (?!\w)  # use negative lookahead since '.' confuses \b test
503    """, re.VERBOSE)
504def _tag_float(source, constant_flag):
505    # Convert floating point constants to single by adding 'f' to the end,
506    # or long double with an 'L' suffix.  OS/X complains if you don't do this.
507    out = FLOAT_RE.sub(r'\g<0>%s'%constant_flag, source)
508    #print("in",repr(source),"out",repr(out), constant_flag)
509    return out
510
511def test_tag_float():
512    """check that floating point constants are properly identified and tagged with 'f'"""
513
514    cases = """
515ZP  : 0.
516ZPF : 0.0,0.01,0.1
517Z  E: 0e+001
518ZP E: 0.E0
519ZPFE: 0.13e-031
520NP  : 1., 12.
521NPF : 1.0001, 1.1, 1.0
522N  E: 1e0, 37E-080
523NP E: 1.e0, 37.E-080
524NPFE: 845.017e+22
525 PF : .1, .0, .0100
526 PFE: .6e+9, .82E-004
527# isolated cases
5280.
5291e0
5300.13e-013
531# untouched
532struct3.e3, 03.05.67, 37
533# expressions
5343.75+-1.6e-7-27+13.2
535a3.e2 - 0.
5364*atan(1)
5374.*atan(1.)
538"""
539
540    output = """
541ZP  : 0.f
542ZPF : 0.0f,0.01f,0.1f
543Z  E: 0e+001f
544ZP E: 0.E0f
545ZPFE: 0.13e-031f
546NP  : 1.f, 12.f
547NPF : 1.0001f, 1.1f, 1.0f
548N  E: 1e0f, 37E-080f
549NP E: 1.e0f, 37.E-080f
550NPFE: 845.017e+22f
551 PF : .1f, .0f, .0100f
552 PFE: .6e+9f, .82E-004f
553# isolated cases
5540.f
5551e0f
5560.13e-013f
557# untouched
558struct3.e3, 03.05.67, 37
559# expressions
5603.75f+-1.6e-7f-27+13.2f
561a3.e2 - 0.f
5624*atan(1)
5634.f*atan(1.f)
564"""
565
566    for case_in, case_out in zip(cases.split('\n'), output.split('\n')):
567        out = _tag_float(case_in, 'f')
568        assert case_out == out, "%r => %r"%(case_in, out)
569
570
571def kernel_name(model_info, variant):
572    # type: (ModelInfo, str) -> str
573    """
574    Name of the exported kernel symbol.
575
576    *variant* is "Iq", "Iqxy" or "Imagnetic".
577    """
578    return model_info.name + "_" + variant
579
580
581def indent(s, depth):
582    # type: (str, int) -> str
583    """
584    Indent a string of text with *depth* additional spaces on each line.
585    """
586    spaces = " "*depth
587    sep = "\n" + spaces
588    return spaces + sep.join(s.split("\n"))
589
590
591_template_cache = {}  # type: Dict[str, Tuple[int, str, str]]
592def load_template(filename):
593    # type: (str) -> str
594    path = joinpath(DATA_PATH, filename)
595    mtime = getmtime(path)
596    if filename not in _template_cache or mtime > _template_cache[filename][0]:
597        with open(path) as fid:
598            _template_cache[filename] = (mtime, fid.read(), path)
599    return _template_cache[filename][1], path
600
601
602_FN_TEMPLATE = """\
603double %(name)s(%(pars)s);
604double %(name)s(%(pars)s) {
605#line %(line)d "%(filename)s"
606    %(body)s
607}
608
609"""
610def _gen_fn(name, pars, body, filename, line):
611    # type: (str, List[Parameter], str, str, int) -> str
612    """
613    Generate a function given pars and body.
614
615    Returns the following string::
616
617         double fn(double a, double b, ...);
618         double fn(double a, double b, ...) {
619             ....
620         }
621    """
622    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
623    return _FN_TEMPLATE % {
624        'name': name, 'pars': par_decl, 'body': body,
625        'filename': filename.replace('\\', '\\\\'), 'line': line,
626    }
627
628
629def _call_pars(prefix, pars):
630    # type: (str, List[Parameter]) -> List[str]
631    """
632    Return a list of *prefix+parameter* from parameter items.
633
634    *prefix* should be "v." if v is a struct.
635    """
636    return [p.as_call_reference(prefix) for p in pars]
637
638
639# type in IQXY pattern could be single, float, double, long double, ...
640_IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)",
641                           flags=re.MULTILINE)
642def _have_Iqxy(sources):
643    # type: (List[str]) -> bool
644    """
645    Return true if any file defines Iqxy.
646
647    Note this is not a C parser, and so can be easily confused by
648    non-standard syntax.  Also, it will incorrectly identify the following
649    as having Iqxy::
650
651        /*
652        double Iqxy(qx, qy, ...) { ... fill this in later ... }
653        */
654
655    If you want to comment out an Iqxy function, use // on the front of the
656    line instead.
657    """
658    for _path, code in sources:
659        if _IQXY_PATTERN.search(code):
660            return True
661    return False
662
663
664def _add_source(source, code, path):
665    """
666    Add a file to the list of source code chunks, tagged with path and line.
667    """
668    path = path.replace('\\', '\\\\')
669    source.append('#line 1 "%s"' % path)
670    source.append(code)
671
672def make_source(model_info):
673    # type: (ModelInfo) -> Dict[str, str]
674    """
675    Generate the OpenCL/ctypes kernel from the module info.
676
677    Uses source files found in the given search path.  Returns None if this
678    is a pure python model, with no C source components.
679    """
680    if callable(model_info.Iq):
681        raise ValueError("can't compile python model")
682        #return None
683
684    # TODO: need something other than volume to indicate dispersion parameters
685    # No volume normalization despite having a volume parameter.
686    # Thickness is labelled a volume in order to trigger polydispersity.
687    # May want a separate dispersion flag, or perhaps a separate category for
688    # disperse, but not volume.  Volume parameters also use relative values
689    # for the distribution rather than the absolute values used by angular
690    # dispersion.  Need to be careful that necessary parameters are available
691    # for computing volume even if we allow non-disperse volume parameters.
692
693    partable = model_info.parameters
694
695    # Load templates and user code
696    kernel_header = load_template('kernel_header.c')
697    kernel_code = load_template('kernel_iq.c')
698    user_code = [(f, open(f).read()) for f in model_sources(model_info)]
699
700    # What kind of 2D model do we need?
701    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str)
702               else 'qac' if not partable.is_asymmetric
703               else 'qabc')
704
705    # Build initial sources
706    source = []
707    _add_source(source, *kernel_header)
708    for path, code in user_code:
709        _add_source(source, code, path)
710
711    # Make parameters for q, qx, qy so that we can use them in declarations
712    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')]
713    # Generate form_volume function, etc. from body only
714    if isinstance(model_info.form_volume, str):
715        pars = partable.form_volume_parameters
716        source.append(_gen_fn('form_volume', pars, model_info.form_volume,
717                              model_info.filename, model_info._form_volume_line))
718    if isinstance(model_info.Iq, str):
719        pars = [q] + partable.iq_parameters
720        source.append(_gen_fn('Iq', pars, model_info.Iq,
721                              model_info.filename, model_info._Iq_line))
722    if isinstance(model_info.Iqxy, str):
723        pars = [qx, qy] + partable.iqxy_parameters
724        source.append(_gen_fn('Iqxy', pars, model_info.Iqxy,
725                              model_info.filename, model_info._Iqxy_line))
726
727    # Define the parameter table
728    lineno = getframeinfo(currentframe()).lineno + 2
729    source.append('#line %d "sasmodels/generate.py"'%lineno)
730    #source.append('introduce breakage in generate to test lineno reporting')
731    source.append("#define PARAMETER_TABLE \\")
732    source.append("\\\n".join(p.as_definition()
733                              for p in partable.kernel_parameters))
734
735    # Define the function calls
736    if partable.form_volume_parameters:
737        refs = _call_pars("_v.", partable.form_volume_parameters)
738        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs))
739    else:
740        # Model doesn't have volume.  We could make the kernel run a little
741        # faster by not using/transferring the volume normalizations, but
742        # the ifdef's reduce readability more than is worthwhile.
743        call_volume = "#define CALL_VOLUME(v) 1.0"
744    source.append(call_volume)
745
746    model_refs = _call_pars("_v.", partable.iq_parameters)
747    pars = ",".join(["_q"] + model_refs)
748    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
749    if xy_mode == 'qabc':
750        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
751        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars
752        clear_iqxy = "#undef CALL_IQ_ABC"
753    elif xy_mode == 'qac':
754        pars = ",".join(["_qa", "_qc"] + model_refs)
755        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars
756        clear_iqxy = "#undef CALL_IQ_AC"
757    else:  # xy_mode == 'qa'
758        pars = ",".join(["_qa"] + model_refs)
759        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
760        clear_iqxy = "#undef CALL_IQ_A"
761
762    magpars = [k-2 for k, p in enumerate(partable.call_parameters)
763               if p.type == 'sld']
764
765    # Fill in definitions for numbers of parameters
766    source.append("#define MAX_PD %s"%partable.max_pd)
767    source.append("#define NUM_PARS %d"%partable.npars)
768    source.append("#define NUM_VALUES %d" % partable.nvalues)
769    source.append("#define NUM_MAGNETIC %d" % partable.nmagnetic)
770    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
771    source.append("#define PROJECTION %d"%PROJECTION)
772
773    # TODO: allow mixed python/opencl kernels?
774
775    ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
776    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
777    result = {
778        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
779        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
780    }
781
782    return result
783
784
785def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name):
786    # type: ([str,str], str, str, str) -> List[str]
787    code = kernel[0]
788    path = kernel[1].replace('\\', '\\\\')
789    iq = [
790        # define the Iq kernel
791        "#define KERNEL_NAME %s_Iq" % name,
792        call_iq,
793        '#line 1 "%s Iq"' % path,
794        code,
795        "#undef CALL_IQ",
796        "#undef KERNEL_NAME",
797        ]
798
799    iqxy = [
800        # define the Iqxy kernel from the same source with different #defines
801        "#define KERNEL_NAME %s_Iqxy" % name,
802        call_iqxy,
803        '#line 1 "%s Iqxy"' % path,
804        code,
805        clear_iqxy,
806        "#undef KERNEL_NAME",
807    ]
808
809    imagnetic = [
810        # define the Imagnetic kernel
811        "#define KERNEL_NAME %s_Imagnetic" % name,
812        "#define MAGNETIC 1",
813        call_iqxy,
814        '#line 1 "%s Imagnetic"' % path,
815        code,
816        clear_iqxy,
817        "#undef MAGNETIC",
818        "#undef KERNEL_NAME",
819    ]
820
821    return iq, iqxy, imagnetic
822
823
824def load_kernel_module(model_name):
825    # type: (str) -> module
826    """
827    Return the kernel module named in *model_name*.
828
829    If the name ends in *.py* then load it as a custom model using
830    :func:`sasmodels.custom.load_custom_kernel_module`, otherwise
831    load it from :mod:`sasmodels.models`.
832    """
833    if model_name.endswith('.py'):
834        kernel_module = load_custom_kernel_module(model_name)
835    else:
836        from sasmodels import models
837        __import__('sasmodels.models.'+model_name)
838        kernel_module = getattr(models, model_name, None)
839    return kernel_module
840
841
842section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z'
843                            % re.escape(string.punctuation))
844def _convert_section_titles_to_boldface(lines):
845    # type: (Sequence[str]) -> Iterator[str]
846    """
847    Do the actual work of identifying and converting section headings.
848    """
849    prior = None
850    for line in lines:
851        if prior is None:
852            prior = line
853        elif section_marker.match(line):
854            if len(line) >= len(prior):
855                yield "".join(("**", prior, "**"))
856                prior = None
857            else:
858                yield prior
859                prior = line
860        else:
861            yield prior
862            prior = line
863    if prior is not None:
864        yield prior
865
866
867def convert_section_titles_to_boldface(s):
868    # type: (str) -> str
869    """
870    Use explicit bold-face rather than section headings so that the table of
871    contents is not polluted with section names from the model documentation.
872
873    Sections are identified as the title line followed by a line of punctuation
874    at least as long as the title line.
875    """
876    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
877
878
879def make_doc(model_info):
880    # type: (ModelInfo) -> str
881    """
882    Return the documentation for the model.
883    """
884    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
885    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
886    docs = model_info.docs if model_info.docs is not None else ""
887    docs = convert_section_titles_to_boldface(docs)
888    pars = make_partable(model_info.parameters.COMMON
889                         + model_info.parameters.kernel_parameters)
890    subst = dict(id=model_info.id.replace('_', '-'),
891                 name=model_info.name,
892                 title=model_info.title,
893                 parameters=pars,
894                 returns=Sq_units if model_info.structure_factor else Iq_units,
895                 docs=docs)
896    return DOC_HEADER % subst
897
898
899# TODO: need a single source for rst_prolog; it is also in doc/rst_prolog
900RST_PROLOG = r"""\
901.. |Ang| unicode:: U+212B
902.. |Ang^-1| replace:: |Ang|\ :sup:`-1`
903.. |Ang^2| replace:: |Ang|\ :sup:`2`
904.. |Ang^-2| replace:: |Ang|\ :sup:`-2`
905.. |1e-6Ang^-2| replace:: 10\ :sup:`-6`\ |Ang|\ :sup:`-2`
906.. |Ang^3| replace:: |Ang|\ :sup:`3`
907.. |Ang^-3| replace:: |Ang|\ :sup:`-3`
908.. |Ang^-4| replace:: |Ang|\ :sup:`-4`
909.. |cm^-1| replace:: cm\ :sup:`-1`
910.. |cm^2| replace:: cm\ :sup:`2`
911.. |cm^-2| replace:: cm\ :sup:`-2`
912.. |cm^3| replace:: cm\ :sup:`3`
913.. |1e15cm^3| replace:: 10\ :sup:`15`\ cm\ :sup:`3`
914.. |cm^-3| replace:: cm\ :sup:`-3`
915.. |sr^-1| replace:: sr\ :sup:`-1`
916
917.. |cdot| unicode:: U+00B7
918.. |deg| unicode:: U+00B0
919.. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3`
920.. |mg/m^2| replace:: mg\ |cdot|\ m\ :sup:`-2`
921.. |fm^2| replace:: fm\ :sup:`2`
922.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1`
923"""
924
925# TODO: make a better fake reference role
926RST_ROLES = """\
927.. role:: ref
928
929.. role:: numref
930
931"""
932
933def make_html(model_info):
934    # type: (ModelInfo) -> str
935    """
936    Convert model docs directly to html.
937    """
938    from . import rst2html
939
940    rst = make_doc(model_info)
941    return rst2html.rst2html("".join((RST_ROLES, RST_PROLOG, rst)))
942
943def view_html(model_name):
944    # type: (str) -> None
945    """
946    Load the model definition and view its help.
947    """
948    from . import modelinfo
949    kernel_module = load_kernel_module(model_name)
950    info = modelinfo.make_model_info(kernel_module)
951    view_html_from_info(info)
952
953def view_html_from_info(info):
954    # type: (ModelInfo) -> None
955    """
956    View the help for a loaded model definition.
957    """
958    from . import rst2html
959    url = "file://"+dirname(info.filename)+"/"
960    rst2html.view_html(make_html(info), url=url)
961
962def demo_time():
963    # type: () -> None
964    """
965    Show how long it takes to process a model.
966    """
967    import datetime
968    from .modelinfo import make_model_info
969    from .models import cylinder
970
971    tic = datetime.datetime.now()
972    make_source(make_model_info(cylinder))
973    toc = (datetime.datetime.now() - tic).total_seconds()
974    print("time: %g"%toc)
975
976
977def main():
978    # type: () -> None
979    """
980    Program which prints the source produced by the model.
981    """
982    from .modelinfo import make_model_info
983
984    if len(sys.argv) <= 1:
985        print("usage: python -m sasmodels.generate modelname")
986    else:
987        name = sys.argv[1]
988        kernel_module = load_kernel_module(name)
989        model_info = make_model_info(kernel_module)
990        source = make_source(model_info)
991        print(source['dll'])
992
993
994if __name__ == "__main__":
995    main()
Note: See TracBrowser for help on using the repository browser.