source: sasmodels/sasmodels/generate.py @ ff10479

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

set projection used in theory to that selected in jitter.py

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