source: sasmodels/sasmodels/generate.py @ d2bb604

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

fix models so all dll tests pass

  • Property mode set to 100644
File size: 31.2 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.
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
125An *model_info* dictionary is constructed from the kernel meta data and
126returned to the caller.
127
128The model evaluator, function call sequence consists of q inputs and the return vector,
129followed by the loop value/weight vector, followed by the values for
130the non-polydisperse parameters, followed by the lengths of the
131polydispersity loops.  To construct the call for 1D models, the
132categories *fixed-1d* and *pd-1d* list the names of the parameters
133of the non-polydisperse and the polydisperse parameters respectively.
134Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models.
135The *pd-rel* category is a set of those parameters which give
136polydispersitiy as a portion of the value (so a 10% length dispersity
137would use a polydispersity value of 0.1) rather than absolute
138dispersity such as an angle plus or minus 15 degrees.
139
140The *volume* category lists the volume parameters in order for calls
141to volume within the kernel (used for volume normalization) and for
142calls to ER and VR for effective radius and volume ratio respectively.
143
144The *orientation* and *magnetic* categories list the orientation and
145magnetic parameters.  These are used by the sasview interface.  The
146blank category is for parameters such as scale which don't have any
147other marking.
148
149The doc string at the start of the kernel module will be used to
150construct the model documentation web pages.  Embedded figures should
151appear in the subdirectory "img" beside the model definition, and tagged
152with the kernel module name to avoid collision with other models.  Some
153file systems are case-sensitive, so only use lower case characters for
154file names and extensions.
155
156
157The function :func:`make` loads the metadata from the module and returns
158the kernel source.  The function :func:`make_doc` extracts the doc string
159and adds the parameter table to the top.  The function :func:`model_sources`
160returns a list of files required by the model.
161
162Code follows the C99 standard with the following extensions and conditions::
163
164    M_PI_180 = pi/180
165    M_4PI_3 = 4pi/3
166    square(x) = x*x
167    cube(x) = x*x*x
168    sinc(x) = sin(x)/x, with sin(0)/0 -> 1
169    all double precision constants must include the decimal point
170    all double declarations may be converted to half, float, or long double
171    FLOAT_SIZE is the number of bytes in the converted variables
172"""
173from __future__ import print_function
174
175#TODO: determine which functions are useful outside of generate
176#__all__ = ["model_info", "make_doc", "make_source", "convert_type"]
177
178from os.path import abspath, dirname, join as joinpath, exists, basename, \
179    splitext, getmtime
180import re
181import string
182import warnings
183
184import numpy as np
185
186from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo
187from .custom import load_custom_kernel_module
188
189TEMPLATE_ROOT = dirname(__file__)
190
191F16 = np.dtype('float16')
192F32 = np.dtype('float32')
193F64 = np.dtype('float64')
194try:  # CRUFT: older numpy does not support float128
195    F128 = np.dtype('float128')
196except TypeError:
197    F128 = None
198
199# Conversion from units defined in the parameter table for each model
200# to units displayed in the sphinx documentation.
201RST_UNITS = {
202    "Ang": "|Ang|",
203    "1/Ang": "|Ang^-1|",
204    "1/Ang^2": "|Ang^-2|",
205    "1e-6/Ang^2": "|1e-6Ang^-2|",
206    "degrees": "degree",
207    "1/cm": "|cm^-1|",
208    "Ang/cm": "|Ang*cm^-1|",
209    "g/cm3": "|g/cm^3|",
210    "mg/m2": "|mg/m^2|",
211    "": "None",
212    }
213
214# Headers for the parameters tables in th sphinx documentation
215PARTABLE_HEADERS = [
216    "Parameter",
217    "Description",
218    "Units",
219    "Default value",
220    ]
221
222# Minimum width for a default value (this is shorter than the column header
223# width, so will be ignored).
224PARTABLE_VALUE_WIDTH = 10
225
226# Documentation header for the module, giving the model name, its short
227# description and its parameter table.  The remainder of the doc comes
228# from the module docstring.
229DOC_HEADER = """.. _%(id)s:
230
231%(name)s
232=======================================================
233
234%(title)s
235
236%(parameters)s
237
238%(returns)s
239
240%(docs)s
241"""
242
243def format_units(units):
244    """
245    Convert units into ReStructured Text format.
246    """
247    return "string" if isinstance(units, list) else RST_UNITS.get(units, units)
248
249def make_partable(pars):
250    """
251    Generate the parameter table to include in the sphinx documentation.
252    """
253    column_widths = [
254        max(len(p.name) for p in pars),
255        max(len(p.description) for p in pars),
256        max(len(format_units(p.units)) for p in pars),
257        PARTABLE_VALUE_WIDTH,
258        ]
259    column_widths = [max(w, len(h))
260                     for w, h in zip(column_widths, PARTABLE_HEADERS)]
261
262    sep = " ".join("="*w for w in column_widths)
263    lines = [
264        sep,
265        " ".join("%-*s" % (w, h)
266                 for w, h in zip(column_widths, PARTABLE_HEADERS)),
267        sep,
268        ]
269    for p in pars:
270        lines.append(" ".join([
271            "%-*s" % (column_widths[0], p.name),
272            "%-*s" % (column_widths[1], p.description),
273            "%-*s" % (column_widths[2], format_units(p.units)),
274            "%*g" % (column_widths[3], p.default),
275            ]))
276    lines.append(sep)
277    return "\n".join(lines)
278
279def _search(search_path, filename):
280    """
281    Find *filename* in *search_path*.
282
283    Raises ValueError if file does not exist.
284    """
285    for path in search_path:
286        target = joinpath(path, filename)
287        if exists(target):
288            return target
289    raise ValueError("%r not found in %s" % (filename, search_path))
290
291
292def model_sources(model_info):
293    """
294    Return a list of the sources file paths for the module.
295    """
296    search_path = [dirname(model_info['filename']),
297                   abspath(joinpath(dirname(__file__), 'models'))]
298    return [_search(search_path, f) for f in model_info['source']]
299
300def timestamp(model_info):
301    """
302    Return a timestamp for the model corresponding to the most recently
303    changed file or dependency.
304    """
305    source_files = (model_sources(model_info)
306                    + model_templates()
307                    + [model_info['filename']])
308    newest = max(getmtime(f) for f in source_files)
309    return newest
310
311def convert_type(source, dtype):
312    """
313    Convert code from double precision to the desired type.
314
315    Floating point constants are tagged with 'f' for single precision or 'L'
316    for long double precision.
317    """
318    if dtype == F16:
319        fbytes = 2
320        source = _convert_type(source, "half", "f")
321    elif dtype == F32:
322        fbytes = 4
323        source = _convert_type(source, "float", "f")
324    elif dtype == F64:
325        fbytes = 8
326        # no need to convert if it is already double
327    elif dtype == F128:
328        fbytes = 16
329        source = _convert_type(source, "long double", "L")
330    else:
331        raise ValueError("Unexpected dtype in source conversion: %s"%dtype)
332    return ("#define FLOAT_SIZE %d\n"%fbytes)+source
333
334
335def _convert_type(source, type_name, constant_flag):
336    """
337    Replace 'double' with *type_name* in *source*, tagging floating point
338    constants with *constant_flag*.
339    """
340    # Convert double keyword to float/long double/half.
341    # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16.
342    # Assume complex numbers are represented as cdouble which is typedef'd
343    # to double2.
344    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
345                    r'\1%s\2'%type_name, source)
346    # Convert floating point constants to single by adding 'f' to the end,
347    # or long double with an 'L' suffix.  OS/X complains if you don't do this.
348    source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?',
349                    r'\g<0>%s'%constant_flag, source)
350    return source
351
352
353def kernel_name(model_info, is_2d):
354    """
355    Name of the exported kernel symbol.
356    """
357    return model_info['name'] + "_" + ("Iqxy" if is_2d else "Iq")
358
359
360def indent(s, depth):
361    """
362    Indent a string of text with *depth* additional spaces on each line.
363    """
364    spaces = " "*depth
365    sep = "\n" + spaces
366    return spaces + sep.join(s.split("\n"))
367
368
369_template_cache = {}
370def load_template(filename):
371    path = joinpath(TEMPLATE_ROOT, filename)
372    mtime = getmtime(path)
373    if filename not in _template_cache or mtime > _template_cache[filename][0]:
374        with open(path) as fid:
375            _template_cache[filename] = (mtime, fid.read(), path)
376    return _template_cache[filename][1]
377
378def model_templates():
379    # TODO: fails DRY; templates are listed in two places.
380    # should instead have model_info contain a list of paths
381    return [joinpath(TEMPLATE_ROOT, filename)
382            for filename in ('kernel_header.c', 'kernel_iq.c')]
383
384
385_FN_TEMPLATE = """\
386double %(name)s(%(pars)s);
387double %(name)s(%(pars)s) {
388    %(body)s
389}
390
391
392"""
393
394def _gen_fn(name, pars, body):
395    """
396    Generate a function given pars and body.
397
398    Returns the following string::
399
400         double fn(double a, double b, ...);
401         double fn(double a, double b, ...) {
402             ....
403         }
404    """
405    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void'
406    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl}
407
408def _call_pars(prefix, pars):
409    """
410    Return a list of *prefix.parameter* from parameter items.
411    """
412    return [p.as_call_reference(prefix) for p in pars]
413
414_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)",
415                           flags=re.MULTILINE)
416def _have_Iqxy(sources):
417    """
418    Return true if any file defines Iqxy.
419
420    Note this is not a C parser, and so can be easily confused by
421    non-standard syntax.  Also, it will incorrectly identify the following
422    as having Iqxy::
423
424        /*
425        double Iqxy(qx, qy, ...) { ... fill this in later ... }
426        */
427
428    If you want to comment out an Iqxy function, use // on the front of the
429    line instead.
430    """
431    for code in sources:
432        if _IQXY_PATTERN.search(code):
433            return True
434    else:
435        return False
436
437def make_source(model_info):
438    """
439    Generate the OpenCL/ctypes kernel from the module info.
440
441    Uses source files found in the given search path.
442    """
443    if callable(model_info['Iq']):
444        return None
445
446    # TODO: need something other than volume to indicate dispersion parameters
447    # No volume normalization despite having a volume parameter.
448    # Thickness is labelled a volume in order to trigger polydispersity.
449    # May want a separate dispersion flag, or perhaps a separate category for
450    # disperse, but not volume.  Volume parameters also use relative values
451    # for the distribution rather than the absolute values used by angular
452    # dispersion.  Need to be careful that necessary parameters are available
453    # for computing volume even if we allow non-disperse volume parameters.
454
455    partable = model_info['parameters']
456
457    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume.
458    # Note that scale and volume are not possible types.
459
460    # Load templates and user code
461    kernel_header = load_template('kernel_header.c')
462    kernel_code = load_template('kernel_iq.c')
463    user_code = [open(f).read() for f in model_sources(model_info)]
464
465    # Build initial sources
466    source = [kernel_header] + user_code
467
468    # Make parameters for q, qx, qy so that we can use them in declarations
469    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')]
470    # Generate form_volume function, etc. from body only
471    if model_info['form_volume'] is not None:
472        pars = partable.form_volume_parameters
473        source.append(_gen_fn('form_volume', pars, model_info['form_volume']))
474    if model_info['Iq'] is not None:
475        pars = [q] + partable.iq_parameters
476        source.append(_gen_fn('Iq', pars, model_info['Iq']))
477    if model_info['Iqxy'] is not None:
478        pars = [qx, qy] + partable.iqxy_parameters
479        source.append(_gen_fn('Iqxy', pars, model_info['Iqxy']))
480
481    # Define the parameter table
482    source.append("#define PARAMETER_TABLE \\")
483    source.append("\\\n".join(p.as_definition()
484                              for p in partable.kernel_parameters))
485
486    # Define the function calls
487    if partable.form_volume_parameters:
488        refs = _call_pars("_v.", partable.form_volume_parameters)
489        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs))
490    else:
491        # Model doesn't have volume.  We could make the kernel run a little
492        # faster by not using/transferring the volume normalizations, but
493        # the ifdef's reduce readability more than is worthwhile.
494        call_volume = "#define CALL_VOLUME(v) 1.0"
495    source.append(call_volume)
496
497    refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters)
498    call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs))
499    if _have_Iqxy(user_code):
500        # Call 2D model
501        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters)
502        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs))
503    else:
504        # Call 1D model with sqrt(qx^2 + qy^2)
505        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))")
506        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters)
507        pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:]
508        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt))
509
510    # Fill in definitions for numbers of parameters
511    source.append("#define MAX_PD %s"%partable.max_pd)
512    source.append("#define NPARS %d"%partable.npars)
513
514    # TODO: allow mixed python/opencl kernels?
515
516    # define the Iq kernel
517    source.append("#define KERNEL_NAME %s_Iq"%model_info['name'])
518    source.append(call_iq)
519    source.append(kernel_code)
520    source.append("#undef CALL_IQ")
521    source.append("#undef KERNEL_NAME")
522
523    # define the Iqxy kernel from the same source with different #defines
524    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name'])
525    source.append(call_iqxy)
526    source.append(kernel_code)
527    source.append("#undef CALL_IQ")
528    source.append("#undef KERNEL_NAME")
529
530    return '\n'.join(source)
531
532class CoordinationDetails(object):
533    def __init__(self, model_info):
534        parameters = model_info['parameters']
535        max_pd = parameters.max_pd
536        npars = parameters.npars
537        par_offset = 4*max_pd
538        self.details = np.zeros(par_offset + 3*npars + 4, 'i4')
539
540        # generate views on different parts of the array
541        self._pd_par     = self.details[0*max_pd:1*max_pd]
542        self._pd_length  = self.details[1*max_pd:2*max_pd]
543        self._pd_offset  = self.details[2*max_pd:3*max_pd]
544        self._pd_stride  = self.details[3*max_pd:4*max_pd]
545        self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars]
546        self._par_coord  = self.details[par_offset+1*npars:par_offset+2*npars]
547        self._pd_coord   = self.details[par_offset+2*npars:par_offset+3*npars]
548
549        # theta_par is fixed
550        self.details[-1] = parameters.theta_offset
551
552    @property
553    def ctypes(self): return self.details.ctypes
554    @property
555    def pd_par(self): return self._pd_par
556    @property
557    def pd_length(self): return self._pd_length
558    @property
559    def pd_offset(self): return self._pd_offset
560    @property
561    def pd_stride(self): return self._pd_stride
562    @property
563    def pd_coord(self): return self._pd_coord
564    @property
565    def par_coord(self): return self._par_coord
566    @property
567    def par_offset(self): return self._par_offset
568    @property
569    def num_coord(self): return self.details[-2]
570    @num_coord.setter
571    def num_coord(self, v): self.details[-2] = v
572    @property
573    def total_pd(self): return self.details[-3]
574    @total_pd.setter
575    def total_pd(self, v): self.details[-3] = v
576    @property
577    def num_active(self): return self.details[-4]
578    @num_active.setter
579    def num_active(self, v): self.details[-4] = v
580
581    def show(self):
582        print("total_pd", self.total_pd)
583        print("num_active", self.num_active)
584        print("pd_par", self.pd_par)
585        print("pd_length", self.pd_length)
586        print("pd_offset", self.pd_offset)
587        print("pd_stride", self.pd_stride)
588        print("par_offsets", self.par_offset)
589        print("num_coord", self.num_coord)
590        print("par_coord", self.par_coord)
591        print("pd_coord", self.pd_coord)
592        print("theta par", self.details[-1])
593
594def mono_details(model_info):
595    details = CoordinationDetails(model_info)
596    # The zero defaults for monodisperse systems are mostly fine
597    details.par_offset[:] = np.arange(2, len(details.par_offset)+2)
598    return details
599
600def poly_details(model_info, weights):
601    #print("weights",weights)
602    weights = weights[2:] # Skip scale and background
603
604    # Decreasing list of polydispersity lengths
605    # Note: the reversing view, x[::-1], does not require a copy
606    pd_length = np.array([len(w) for w in weights])
607    num_active = np.sum(pd_length>1)
608    if num_active > model_info['parameters'].max_pd:
609        raise ValueError("Too many polydisperse parameters")
610
611    pd_offset = np.cumsum(np.hstack((0, pd_length)))
612    idx = np.argsort(pd_length)[::-1][:num_active]
613    par_length = np.array([max(len(w),1) for w in weights])
614    pd_stride = np.cumprod(np.hstack((1, par_length[idx])))
615    par_offsets = np.cumsum(np.hstack((2, par_length)))
616
617    details = CoordinationDetails(model_info)
618    details.pd_par[:num_active] = idx
619    details.pd_length[:num_active] = pd_length[idx]
620    details.pd_offset[:num_active] = pd_offset[idx]
621    details.pd_stride[:num_active] = pd_stride[:-1]
622    details.par_offset[:] = par_offsets[:-1]
623    details.total_pd = pd_stride[-1]
624    details.num_active = num_active
625    # Without constraints coordinated parameters are just the pd parameters
626    details.par_coord[:num_active] = idx
627    details.pd_coord[:num_active] = 2**np.arange(num_active)
628    details.num_coord = num_active
629    #details.show()
630    return details
631
632def constrained_poly_details(model_info, weights, constraints):
633    # Need to find the independently varying pars and sort them
634    # Need to build a coordination list for the dependent variables
635    # Need to generate a constraints function which takes values
636    # and weights, returning par blocks
637    raise NotImplementedError("Can't handle constraints yet")
638
639
640def create_default_functions(model_info):
641    """
642    Autogenerate missing functions, such as Iqxy from Iq.
643
644    This only works for Iqxy when Iq is written in python. :func:`make_source`
645    performs a similar role for Iq written in C.
646    """
647    if callable(model_info['Iq']) and model_info['Iqxy'] is None:
648        partable = model_info['parameters']
649        if partable.has_2d:
650            raise ValueError("Iqxy model is missing")
651        Iq = model_info['Iq']
652        def Iqxy(qx, qy, **kw):
653            return Iq(np.sqrt(qx**2 + qy**2), **kw)
654        model_info['Iqxy'] = Iqxy
655
656
657def load_kernel_module(model_name):
658    if model_name.endswith('.py'):
659        kernel_module = load_custom_kernel_module(model_name)
660    else:
661        from sasmodels import models
662        __import__('sasmodels.models.'+model_name)
663        kernel_module = getattr(models, model_name, None)
664    return kernel_module
665
666
667def make_model_info(kernel_module):
668    """
669    Interpret the model definition file, categorizing the parameters.
670
671    The module can be loaded with a normal python import statement if you
672    know which module you need, or with __import__('sasmodels.model.'+name)
673    if the name is in a string.
674
675    The *model_info* structure contains the following fields:
676
677    * *id* is the id of the kernel
678    * *name* is the display name of the kernel
679    * *filename* is the full path to the module defining the file (if any)
680    * *title* is a short description of the kernel
681    * *description* is a long description of the kernel (this doesn't seem
682      very useful since the Help button on the model page brings you directly
683      to the documentation page)
684    * *docs* is the docstring from the module.  Use :func:`make_doc` to
685    * *category* specifies the model location in the docs
686    * *parameters* is the model parameter table
687    * *single* is True if the model allows single precision
688    * *structure_factor* is True if the model is useable in a product
689    * *variant_info* contains the information required to select between
690      model variants (e.g., the list of cases) or is None if there are no
691      model variants
692    * *par_type* categorizes the model parameters. See
693      :func:`categorize_parameters` for details.
694    * *demo* contains the *{parameter: value}* map used in compare (and maybe
695      for the demo plot, if plots aren't set up to use the default values).
696      If *demo* is not given in the file, then the default values will be used.
697    * *tests* is a set of tests that must pass
698    * *source* is the list of library files to include in the C model build
699    * *Iq*, *Iqxy*, *form_volume*, *ER*, *VR* and *sesans* are python functions
700      implementing the kernel for the module, or None if they are not
701      defined in python
702    * *composition* is None if the model is independent, otherwise it is a
703      tuple with composition type ('product' or 'mixture') and a list of
704      *model_info* blocks for the composition objects.  This allows us to
705      build complete product and mixture models from just the info.
706    """
707    # TODO: maybe turn model_info into a class ModelDefinition
708    #print("make parameter table", kernel_module.parameters)
709    parameters = make_parameter_table(kernel_module.parameters)
710    filename = abspath(kernel_module.__file__)
711    kernel_id = splitext(basename(filename))[0]
712    name = getattr(kernel_module, 'name', None)
713    if name is None:
714        name = " ".join(w.capitalize() for w in kernel_id.split('_'))
715    model_info = dict(
716        id=kernel_id,  # string used to load the kernel
717        filename=abspath(kernel_module.__file__),
718        name=name,
719        title=getattr(kernel_module, 'title', name+" model"),
720        description=getattr(kernel_module, 'description', 'no description'),
721        parameters=parameters,
722        composition=None,
723        docs=kernel_module.__doc__,
724        category=getattr(kernel_module, 'category', None),
725        single=getattr(kernel_module, 'single', True),
726        structure_factor=getattr(kernel_module, 'structure_factor', False),
727        profile_axes=getattr(kernel_module, 'profile_axes', ['x','y']),
728        variant_info=getattr(kernel_module, 'invariant_info', None),
729        demo=getattr(kernel_module, 'demo', None),
730        source=getattr(kernel_module, 'source', []),
731        tests=getattr(kernel_module, 'tests', []),
732        )
733    set_demo(model_info, getattr(kernel_module, 'demo', None))
734    # Check for optional functions
735    functions = "ER VR form_volume Iq Iqxy profile sesans".split()
736    model_info.update((k, getattr(kernel_module, k, None)) for k in functions)
737    create_default_functions(model_info)
738    # Precalculate the monodisperse parameters
739    # TODO: make this a lazy evaluator
740    # make_model_info is called for every model on sasview startup
741    model_info['mono_details'] = mono_details(model_info)
742    return model_info
743
744section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z'
745                            %re.escape(string.punctuation))
746def _convert_section_titles_to_boldface(lines):
747    """
748    Do the actual work of identifying and converting section headings.
749    """
750    prior = None
751    for line in lines:
752        if prior is None:
753            prior = line
754        elif section_marker.match(line):
755            if len(line) >= len(prior):
756                yield "".join(("**", prior, "**"))
757                prior = None
758            else:
759                yield prior
760                prior = line
761        else:
762            yield prior
763            prior = line
764    if prior is not None:
765        yield prior
766
767def convert_section_titles_to_boldface(s):
768    """
769    Use explicit bold-face rather than section headings so that the table of
770    contents is not polluted with section names from the model documentation.
771
772    Sections are identified as the title line followed by a line of punctuation
773    at least as long as the title line.
774    """
775    return "\n".join(_convert_section_titles_to_boldface(s.split('\n')))
776
777def make_doc(model_info):
778    """
779    Return the documentation for the model.
780    """
781    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale."
782    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$."
783    docs = convert_section_titles_to_boldface(model_info['docs'])
784    subst = dict(id=model_info['id'].replace('_', '-'),
785                 name=model_info['name'],
786                 title=model_info['title'],
787                 parameters=make_partable(model_info['parameters']),
788                 returns=Sq_units if model_info['structure_factor'] else Iq_units,
789                 docs=docs)
790    return DOC_HEADER % subst
791
792
793def demo_time():
794    """
795    Show how long it takes to process a model.
796    """
797    from .models import cylinder
798    import datetime
799    tic = datetime.datetime.now()
800    make_source(make_model_info(cylinder))
801    toc = (datetime.datetime.now() - tic).total_seconds()
802    print("time: %g"%toc)
803
804def main():
805    """
806    Program which prints the source produced by the model.
807    """
808    import sys
809    if len(sys.argv) <= 1:
810        print("usage: python -m sasmodels.generate modelname")
811    else:
812        name = sys.argv[1]
813        kernel_module = load_kernel_module(name)
814        model_info = make_model_info(kernel_module)
815        source = make_source(model_info)
816        print(source)
817
818if __name__ == "__main__":
819    main()
Note: See TracBrowser for help on using the repository browser.