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

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

fix py3 support for opencl kernels

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