source: sasmodels/sasmodels/generate.py @ 1e7b0db0

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 1e7b0db0 was 1e7b0db0, checked in by wojciech, 6 years ago

sinc tranfered to sas_sinx_x

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