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

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 2e49f9e was 2e49f9e, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

improve plugin portability: OSX wants 1e8⇒1e8f and acos(0)⇒ acos(0.f)

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