source: sasmodels/sasmodels/modelinfo.py @ b2f0e5d

ticket-1257-vesicle-product
Last change on this file since b2f0e5d was b2f0e5d, checked in by Paul Kienzle <pkienzle@…>, 7 months ago

Fix handling of volfraction in vesicle@hardsphere. Refs 1257.

  • Property mode set to 100644
File size: 48.0 KB
Line 
1"""
2Model Info and Parameter Tables
3===============================
4
5Defines :class:`ModelInfo` and :class:`ParameterTable` and the routines for
6manipulating them.  In particular, :func:`make_model_info` converts a kernel
7module into the model info block as seen by the rest of the sasmodels library.
8"""
9from __future__ import print_function
10
11from copy import copy
12from os.path import abspath, basename, splitext
13import inspect
14
15import numpy as np  # type: ignore
16
17# Optional typing
18# pylint: disable=unused-import
19try:
20    from typing import Tuple, List, Union, Dict, Optional, Any, Callable, Sequence, Set
21    from types import ModuleType
22except ImportError:
23    pass
24else:
25    Limits = Tuple[float, float]
26    #LimitsOrChoice = Union[Limits, Tuple[Sequence[str]]]
27    ParameterDef = Tuple[str, str, float, Limits, str, str]
28    ParameterSetUser = Dict[str, Union[float, List[float]]]
29    ParameterSet = Dict[str, float]
30    TestInput = Union[str, float, List[float], Tuple[float, float], List[Tuple[float, float]]]
31    TestValue = Union[float, List[float]]
32    TestCondition = Tuple[ParameterSetUser, TestInput, TestValue]
33# pylint: enable=unused-import
34
35# If MAX_PD changes, need to change the loop macros in kernel_iq.c
36MAX_PD = 5 #: Maximum number of simultaneously polydisperse parameters
37
38# assumptions about common parameters exist throughout the code, such as:
39# (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them
40# (2) kernel drivers assume scale is par[0] and background is par[1]
41# (3) mixture models drop the background on components and replace the scale
42#     with a scale that varies from [-inf, inf]
43# (4) product models drop the background and reassign scale
44# and maybe other places.
45# Note that scale and background cannot be coordinated parameters whose value
46# depends on the some polydisperse parameter with the current implementation
47DEFAULT_BACKGROUND = 1e-3
48COMMON_PARAMETERS = [
49    ("scale", "", 1, (0.0, np.inf), "", "Source intensity"),
50    ("background", "1/cm", DEFAULT_BACKGROUND, (-np.inf, np.inf), "", "Source background"),
51]
52assert (len(COMMON_PARAMETERS) == 2
53        and COMMON_PARAMETERS[0][0] == "scale"
54        and COMMON_PARAMETERS[1][0] == "background"), "don't change common parameters"
55
56
57def make_parameter_table(pars):
58    # type: (List[ParameterDef]) -> ParameterTable
59    """
60    Construct a parameter table from a list of parameter definitions.
61
62    This is used by the module processor to convert the parameter block into
63    the parameter table seen in the :class:`ModelInfo` for the module.
64    """
65    processed = []
66    for p in pars:
67        if not isinstance(p, (list, tuple)) or len(p) != 6:
68            raise ValueError("Parameter should be [name, units, default, limits, type, desc], but got %r"
69                             %str(p))
70        processed.append(parse_parameter(*p))
71    partable = ParameterTable(processed)
72    partable.check_angles(strict=True)
73    return partable
74
75def parse_parameter(name, units='', default=np.NaN,
76                    user_limits=None, ptype='', description=''):
77    # type: (str, str, float, Sequence[Any], str, str) -> Parameter
78    """
79    Parse an individual parameter from the parameter definition block.
80
81    This does type and value checking on the definition, leading
82    to early failure in the model loading process and easier debugging.
83    """
84    # Parameter is a user facing class.  Do robust type checking.
85    if not isstr(name):
86        raise ValueError("expected string for parameter name %r"%name)
87    if not isstr(units):
88        raise ValueError("expected units to be a string for %s"%name)
89
90    # Process limits as [float, float] or [[str, str, ...]]
91    choices = []  # type: List[str]
92    if user_limits is None:
93        limits = (-np.inf, np.inf)
94    elif not isinstance(user_limits, (tuple, list)):
95        raise ValueError("invalid limits for %s"%name)
96    else:
97        # if limits is [[str,...]], then this is a choice list field,
98        # and limits are 1 to length of string list
99        if isinstance(user_limits[0], (tuple, list)):
100            choices = user_limits[0]
101            limits = (0., len(choices)-1.)
102            if not all(isstr(k) for k in choices):
103                raise ValueError("choices must be strings for %s"%name)
104        else:
105            try:
106                low, high = user_limits
107                limits = (float(low), float(high))
108            except Exception:
109                raise ValueError("invalid limits for %s: %r"%(name, user_limits))
110            if low >= high:
111                raise ValueError("require lower limit < upper limit")
112
113    # Process default value as float, making sure it is in range
114    if not isinstance(default, (int, float)):
115        raise ValueError("expected default %r to be a number for %s"
116                         % (default, name))
117    if default < limits[0] or default > limits[1]:
118        raise ValueError("default value %r not in range for %s"
119                         % (default, name))
120
121    # Check for valid parameter type
122    if ptype not in ("volume", "orientation", "sld", "magnetic", ""):
123        raise ValueError("unexpected type %r for %s" % (ptype, name))
124
125    # Check for valid parameter description
126    if not isstr(description):
127        raise ValueError("expected description to be a string")
128
129    # Parameter id for name[n] does not include [n]
130    if "[" in name:
131        if not name.endswith(']'):
132            raise ValueError("Expected name[len] for vector parameter %s"%name)
133        pid, ref = name[:-1].split('[', 1)
134        ref = ref.strip()
135    else:
136        pid, ref = name, None
137
138    # automatically identify sld types
139    if ptype == '' and (pid.startswith('sld') or pid.endswith('sld')):
140        ptype = 'sld'
141
142    # Check if using a vector definition, name[k], as the parameter name
143    if ref:
144        if ref == '':
145            raise ValueError("Need to specify vector length for %s"%name)
146        try:
147            length = int(ref)
148            control = None
149        except ValueError:
150            length = None
151            control = ref
152    else:
153        length = 1
154        control = None
155
156    # Build the parameter
157    parameter = Parameter(name=name, units=units, default=default,
158                          limits=limits, ptype=ptype, description=description)
159
160    # TODO: need better control over whether a parameter is polydisperse
161    parameter.polydisperse = ptype in ('orientation', 'volume')
162    parameter.relative_pd = ptype == 'volume'
163    parameter.choices = choices
164    parameter.length = length
165    parameter.length_control = control
166    return parameter
167
168
169def expand_pars(partable, pars):
170    # type: (ParameterTable, ParameterSetUser) ->  ParameterSet
171    """
172    Create demo parameter set from key-value pairs.
173
174    *pars* are the key-value pairs to use for the parameters.  Any
175    parameters not specified in *pars* are set from the *partable* defaults.
176
177    If *pars* references vector fields, such as thickness[n], then support
178    different ways of assigning the demo values, including assigning a
179    specific value (e.g., thickness3=50.0), assigning a new value to all
180    (e.g., thickness=50.0) or assigning values using list notation.
181    """
182    if pars is None:
183        result = partable.defaults
184    else:
185        lookup = dict((p.id, p) for p in partable.kernel_parameters)
186        result = partable.defaults.copy()
187        scalars = dict((name, value) for name, value in pars.items()
188                       if name not in lookup or lookup[name].length == 1)
189        vectors = dict((name, value) for name, value in pars.items()
190                       if name in lookup and lookup[name].length > 1)
191        #print("lookup", lookup)
192        #print("scalars", scalars)
193        #print("vectors", vectors)
194        if vectors:
195            for name, value in vectors.items():
196                if np.isscalar(value):
197                    # support for the form
198                    #    dict(thickness=0, thickness2=50)
199                    for k in range(1, lookup[name].length+1):
200                        key = name+str(k)
201                        if key not in scalars:
202                            scalars[key] = value
203                else:
204                    # supoprt for the form
205                    #    dict(thickness=[20,10,3])
206                    for (k, v) in enumerate(value):
207                        scalars[name+str(k+1)] = v
208        result.update(scalars)
209        #print("expanded", result)
210
211    return result
212
213def prefix_parameter(par, prefix):
214    # type: (Parameter, str) -> Parameter
215    """
216    Return a copy of the parameter with its name prefixed.
217    """
218    new_par = copy(par)
219    new_par.name = prefix + par.name
220    new_par.id = prefix + par.id
221
222def suffix_parameter(par, suffix):
223    # type: (Parameter, str) -> Parameter
224    """
225    Return a copy of the parameter with its name prefixed.
226    """
227    new_par = copy(par)
228    # If name has the form x[n], replace with x_suffix[n]
229    new_par.name = par.id + suffix + par.name[len(par.id):]
230    new_par.id = par.id + suffix
231
232class Parameter(object):
233    """
234    The available kernel parameters are defined as a list, with each parameter
235    defined as a sublist with the following elements:
236
237    *name* is the name that will be displayed to the user.  Names
238    should be lower case, with words separated by underscore.  If
239    acronyms are used, the whole acronym should be upper case. For vector
240    parameters, the name will be followed by *[len]* where *len* is an
241    integer length of the vector, or the name of the parameter which
242    controls the length.  The attribute *id* will be created from name
243    without the length.
244
245    *units* should be one of *degrees* for angles, *Ang* for lengths,
246    *1e-6/Ang^2* for SLDs.
247
248    *default value* will be the initial value for  the model when it
249    is selected, or when an initial value is not otherwise specified.
250
251    *limits = [lb, ub]* are the hard limits on the parameter value, used to
252    limit the polydispersity density function.  In the fit, the parameter limits
253    given to the fit are the limits  on the central value of the parameter.
254    If there is polydispersity, it will evaluate parameter values outside
255    the fit limits, but not outside the hard limits specified in the model.
256    If there are no limits, use +/-inf imported from numpy.
257
258    *type* indicates how the parameter will be used.  "volume" parameters
259    will be used in all functions.  "orientation" parameters are not passed,
260    but will be used to convert from *qx*, *qy* to *qa*, *qb*, *qc* in calls to
261    *Iqxy* and *Imagnetic*.  If *type* is the empty string, the parameter will
262    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters
263    can automatically be promoted to magnetic parameters, each of which
264    will have a magnitude and a direction, which may be different from
265    other sld parameters. The volume parameters are used for calls
266    to form_volume within the kernel (required for volume normalization),
267    to shell_volume (for hollow shapes), and to radius_effective (for
268    structure factor interactions) respectively.
269
270    *description* is a short description of the parameter.  This will
271    be displayed in the parameter table and used as a tool tip for the
272    parameter value in the user interface.
273
274    Additional values can be set after the parameter is created:
275
276    * *length* is the length of the field if it is a vector field
277
278    * *length_control* is the parameter which sets the vector length
279
280    * *is_control* is True if the parameter is a control parameter for a vector
281
282    * *polydisperse* is true if the parameter accepts a polydispersity
283
284    * *relative_pd* is true if that polydispersity is a portion of the
285      value (so a 10% length dipsersity would use a polydispersity value
286      of 0.1) rather than absolute dispersisity (such as an angle plus or
287      minus 15 degrees).
288
289    *choices* is the option names for a drop down list of options, as for
290    example, might be used to set the value of a shape parameter.
291
292    Control parameters are used for variant models such as :ref:`rpa` which
293    have different cases with different parameters, as well as models
294    like :ref:`spherical_sld` with its user defined number of shells.
295    The control parameter should appear in the parameter table along with the
296    parameters it is is controlling.  For variant models, use *[CASES]* in
297    place of the parameter limits Within the parameter definition table,
298    with case names such as::
299
300         CASES = ["diblock copolymer", "triblock copolymer", ...]
301
302    This should give *limits=[[case1, case2, ...]]*, but the model loader
303    translates it to *limits=[0, len(CASES)-1]*, and adds *choices=CASES* to
304    the :class:`Parameter` definition. Note that models can use a list of
305    cases as a parameter without it being a control parameter.  Either way,
306    the parameter is sent to the model evaluator as *float(choice_num)*,
307    where choices are numbered from 0. :meth:`ModelInfo.get_hidden_parameters`
308    will determine which parameers to display.
309
310    The class contructor should not be called directly, but instead the
311    parameter table is built using :func:`make_parameter_table` and
312    :func:`parse_parameter` therein.
313    """
314    def __init__(self, name, units='', default=None, limits=(-np.inf, np.inf),
315                 ptype='', description=''):
316        # type: (str, str, float, Limits, str, str) -> None
317        self.id = name.split('[')[0].strip() # type: str
318        self.name = name                     # type: str
319        self.units = units                   # type: str
320        self.default = default               # type: float
321        self.limits = limits                 # type: Limits
322        self.type = ptype                    # type: str
323        self.description = description       # type: str
324
325        # Length and length_control will be filled in once the complete
326        # parameter table is available.
327        self.length = 1                      # type: int
328        self.length_control = None           # type: Optional[str]
329        self.is_control = False              # type: bool
330
331        # TODO: need better control over whether a parameter is polydisperse
332        self.polydisperse = False            # type: bool
333        self.relative_pd = False             # type: bool
334
335        # choices are also set externally.
336        self.choices = []                    # type: List[str]
337
338    def as_definition(self):
339        # type: () -> str
340        """
341        Declare space for the variable in a parameter structure.
342
343        For example, the parameter thickness with length 3 will
344        return "double thickness[3];", with no spaces before and
345        no new line character afterward.
346        """
347        if self.length == 1:
348            return "double %s;"%self.id
349        else:
350            return "double %s[%d];"%(self.id, self.length)
351
352    def as_function_argument(self):
353        # type: () -> str
354        r"""
355        Declare the variable as a function argument.
356
357        For example, the parameter thickness with length 3 will
358        return "double \*thickness", with no spaces before and
359        no comma afterward.
360        """
361        if self.length == 1:
362            return "double %s"%self.id
363        else:
364            return "double *%s"%self.id
365
366    def as_call_reference(self, prefix=""):
367        # type: (str) -> str
368        """
369        Return *prefix* + parameter name.  For struct references, use "v."
370        as the prefix.
371        """
372        # Note: if the parameter is a struct type, then we will need to use
373        # &prefix+id.  For scalars and vectors we can just use prefix+id.
374        return prefix + self.id
375
376    def __str__(self):
377        # type: () -> str
378        return "<%s>"%self.name
379
380    def __repr__(self):
381        # type: () -> str
382        return "P<%s>"%self.name
383
384
385class ParameterTable(object):
386    """
387    ParameterTable manages the list of available parameters.
388
389    There are a couple of complications which mean that the list of parameters
390    for the kernel differs from the list of parameters that the user sees.
391
392    (1) Common parameters.  Scale and background are implicit to every model,
393    but are not passed to the kernel.
394
395    (2) Vector parameters.  Vector parameters are passed to the kernel as a
396    pointer to an array, e.g., thick[], but they are seen by the user as n
397    separate parameters thick1, thick2, ...
398
399    Therefore, the parameter table is organized by how it is expected to be
400    used. The following information is needed to set up the kernel functions:
401
402    * *kernel_parameters* is the list of parameters in the kernel parameter
403      table, with vector parameter p declared as p[].
404
405    * *iq_parameters* is the list of parameters to the Iq(q, ...) function,
406      with vector parameter p sent as p[].
407
408    * *form_volume_parameters* is the list of parameters to the form_volume(...)
409      function, with vector parameter p sent as p[].
410
411    Problem details, which sets up the polydispersity loops, requires the
412    following:
413
414    * *theta_offset* is the offset of the theta parameter in the kernel parameter
415      table, with vector parameters counted as n individual parameters
416      p1, p2, ..., or offset is -1 if there is no theta parameter.
417
418    * *max_pd* is the maximum number of polydisperse parameters, with vector
419      parameters counted as n individual parameters p1, p2, ...  Note that
420      this number is limited to sasmodels.modelinfo.MAX_PD.
421
422    * *npars* is the total number of parameters to the kernel, with vector
423      parameters counted as n individual parameters p1, p2, ...
424
425    * *common_parameters* is the list of common parameters, with a unique
426      copy for each model so that structure factors can have a default
427      background of 0.0.
428
429    * *call_parameters* is the complete list of parameters to the kernel,
430      including scale and background, with vector parameters recorded as
431      individual parameters p1, p2, ...
432
433    * *active_1d* is the set of names that may be polydisperse for 1d data
434
435    * *active_2d* is the set of names that may be polydisperse for 2d data
436
437    User parameters are the set of parameters visible to the user, including
438    the scale and background parameters that the kernel does not see.  User
439    parameters don't use vector notation, and instead use p1, p2, ...
440    """
441    def __init__(self, parameters):
442        # type: (List[Parameter]) -> None
443
444        # scale and background are implicit parameters
445        # Need them to be unique to each model in case they have different
446        # properties, such as default=0.0 for structure factor backgrounds.
447        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS]
448        self.kernel_parameters = parameters
449        self._set_vector_lengths()
450        self.npars = sum(p.length for p in self.kernel_parameters)
451        self.nmagnetic = sum(p.length for p in self.kernel_parameters
452                             if p.type == 'sld')
453        self.nvalues = 2 + self.npars
454        if self.nmagnetic:
455            self.nvalues += 3 + 3*self.nmagnetic
456        self.call_parameters = self._get_call_parameters()
457        self.defaults = self._get_defaults()
458        #self._name_table= dict((p.id, p) for p in parameters)
459
460
461        # Set the kernel parameters.  Assumes background and scale are the
462        # first two parameters in the parameter list, but these are not sent
463        # to the underlying kernel functions.
464        self.iq_parameters = [p for p in self.kernel_parameters
465                              if p.type not in ('orientation', 'magnetic')]
466        self.orientation_parameters = [p for p in self.kernel_parameters
467                                       if p.type == 'orientation']
468        self.form_volume_parameters = [p for p in self.kernel_parameters
469                                       if p.type == 'volume']
470
471        # Theta offset
472        offset = 0
473        for p in self.kernel_parameters:
474            if p.name == 'theta':
475                self.theta_offset = offset
476                break
477            offset += p.length
478        else:
479            self.theta_offset = -1
480
481        # number of polydisperse parameters
482        num_pd = sum(p.length for p in self.kernel_parameters if p.polydisperse)
483        # Don't use more polydisperse parameters than are available in the model
484        self.max_pd = min(num_pd, MAX_PD)
485
486        # true if has 2D parameters
487        self.has_2d = any(p.type in ('orientation', 'magnetic')
488                          for p in self.kernel_parameters)
489        self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters)
490        self.magnetism_index = [k for k, p in enumerate(self.call_parameters)
491                                if p.id.endswith('_M0')]
492
493        self.pd_1d = set(p.name for p in self.call_parameters
494                         if p.polydisperse and p.type not in ('orientation', 'magnetic'))
495        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse)
496
497        # Final checks
498        self.check_duplicates()
499        self.check_angles()
500
501    def set_zero_background(self):
502        """
503        Set the default background to zero for this model.  This is done for
504        structure factor models.
505        """
506        # type: () -> None
507        # Make sure background is the second common parameter.
508        assert self.common_parameters[1].id == "background"
509        self.common_parameters[1].default = 0.0
510        self.defaults = self._get_defaults()
511
512    def check_angles(self, strict=False):
513        """
514        Check that orientation angles are theta, phi and possibly psi.
515
516        *strict* should be True when checking a parameter table defined
517        in a model file, but False when checking from mixture models, etc.,
518        where the parameters aren't being passed to a calculator directly.
519        """
520        theta = phi = psi = -1
521        for k, p in enumerate(self.kernel_parameters):
522            if p.name == 'theta':
523                theta = k
524                if p.type != 'orientation':
525                    raise TypeError("theta must be an orientation parameter")
526            elif p.name == 'phi':
527                phi = k
528                if p.type != 'orientation':
529                    raise TypeError("phi must be an orientation parameter")
530            elif p.name == 'psi':
531                psi = k
532                if p.type != 'orientation':
533                    raise TypeError("psi must be an orientation parameter")
534            elif p.type == 'orientation' and strict:
535                raise TypeError("only theta, phi and psi can be orientation parameters")
536        if theta >= 0 and phi >= 0:
537            last_par = len(self.kernel_parameters) - 1
538            if phi != theta+1:
539                raise TypeError("phi must follow theta")
540            if psi >= 0 and psi != phi+1:
541                raise TypeError("psi must follow phi")
542            # TODO: Why must theta/phi/psi be at the end?  Consistency only?
543            if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par):
544                if strict:
545                    raise TypeError("orientation parameters must appear at the "
546                                    "end of the parameter table")
547        elif theta >= 0 or phi >= 0 or psi >= 0:
548            raise TypeError("oriented shapes must have both theta and phi and maybe psi")
549
550    def check_duplicates(self):
551        """
552        Check for duplicate parameter names
553        """
554        checked, dups = set(), set()
555        for p in self.call_parameters:
556            if p.id in checked:
557                dups.add(p.id)
558            else:
559                checked.add(p.id)
560        if dups:
561            raise TypeError("Duplicate parameters: {}"
562                            .format(", ".join(sorted(dups))))
563
564    def __getitem__(self, key):
565        # Find the parameter definition
566        for par in self.call_parameters:
567            if par.name == key:
568                return par
569        raise KeyError("unknown parameter %r"%key)
570
571    def __contains__(self, key):
572        for par in self.call_parameters:
573            if par.name == key:
574                return True
575        return False
576
577    def _set_vector_lengths(self):
578        # type: () -> List[str]
579        """
580        Walk the list of kernel parameters, setting the length field of the
581        vector parameters from the upper limit of the reference parameter.
582
583        This needs to be done once the entire parameter table is available
584        since the reference may still be undefined when the parameter is
585        initially created.
586
587        Returns the list of control parameter names.
588
589        Note: This modifies the underlying parameter object.
590        """
591        # Sort out the length of the vector parameters such as thickness[n]
592        for p in self.kernel_parameters:
593            if p.length_control:
594                ref = self._get_ref(p)
595                ref.is_control = True
596                ref.polydisperse = False
597                low, high = ref.limits
598                if int(low) != low or int(high) != high or low < 0 or high > 20:
599                    raise ValueError("expected limits on %s to be within [0, 20]"
600                                     % ref.name)
601                p.length = int(high)
602
603    def _get_ref(self, p):
604        # type: (Parameter) -> Parameter
605        for ref in self.kernel_parameters:
606            if ref.id == p.length_control:
607                return ref
608        raise ValueError("no reference variable %r for %s"
609                         % (p.length_control, p.name))
610
611    def _get_defaults(self):
612        # type: () -> ParameterSet
613        """
614        Get a list of parameter defaults from the parameters.
615
616        Expands vector parameters into parameter id+number.
617        """
618        # Construct default values, including vector defaults
619        defaults = {}
620        for p in self.call_parameters:
621            if p.length == 1:
622                defaults[p.id] = p.default
623            else:
624                for k in range(1, p.length+1):
625                    defaults["%s%d"%(p.id, k)] = p.default
626        return defaults
627
628    def _get_call_parameters(self):
629        # type: () -> List[Parameter]
630        full_list = self.common_parameters[:]
631        for p in self.kernel_parameters:
632            if p.length == 1:
633                full_list.append(p)
634            else:
635                for k in range(1, p.length+1):
636                    pk = Parameter(p.id+str(k), p.units, p.default,
637                                   p.limits, p.type, p.description)
638                    pk.polydisperse = p.polydisperse
639                    pk.relative_pd = p.relative_pd
640                    pk.choices = p.choices
641                    full_list.append(pk)
642
643        # Add the magnetic parameters to the end of the call parameter list.
644        if self.nmagnetic > 0:
645            full_list.extend([
646                Parameter('up_frac_i', '', 0., [0., 1.],
647                          'magnetic', 'fraction of spin up incident'),
648                Parameter('up_frac_f', '', 0., [0., 1.],
649                          'magnetic', 'fraction of spin up final'),
650                Parameter('up_angle', 'degrees', 0., [0., 360.],
651                          'magnetic', 'spin up angle'),
652            ])
653            slds = [p for p in full_list if p.type == 'sld']
654            for p in slds:
655                full_list.extend([
656                    Parameter(p.id+'_M0', '1e-6/Ang^2', 0., [-np.inf, np.inf],
657                              'magnetic', 'magnetic amplitude for '+p.description),
658                    Parameter(p.id+'_mtheta', 'degrees', 0., [-90., 90.],
659                              'magnetic', 'magnetic latitude for '+p.description),
660                    Parameter(p.id+'_mphi', 'degrees', 0., [-180., 180.],
661                              'magnetic', 'magnetic longitude for '+p.description),
662                ])
663        #print("call parameters", full_list)
664        return full_list
665
666    def user_parameters(self, pars, is2d=True):
667        # type: (Dict[str, float], bool) -> List[Parameter]
668        """
669        Return the list of parameters for the given data type.
670
671        Vector parameters are expanded in place.  If multiple parameters
672        share the same vector length, then the parameters will be interleaved
673        in the result.  The control parameters come first.  For example,
674        if the parameter table is ordered as::
675
676            sld_core
677            sld_shell[num_shells]
678            sld_solvent
679            thickness[num_shells]
680            num_shells
681
682        and *pars[num_shells]=2* then the returned list will be::
683
684            num_shells
685            scale
686            background
687            sld_core
688            sld_shell1
689            thickness1
690            sld_shell2
691            thickness2
692            sld_solvent
693
694        Note that shell/thickness pairs are grouped together in the result
695        even though they were not grouped in the incoming table.  The control
696        parameter is always returned first since the GUI will want to set it
697        early, and rerender the table when it is changed.
698
699        Parameters marked as sld will automatically have a set of associated
700        magnetic parameters (p_M0, p_mtheta, p_mphi), as well as polarization
701        information (up_theta, up_frac_i, up_frac_f).
702        """
703        # control parameters go first
704        control = [p for p in self.kernel_parameters if p.is_control]
705
706        # Gather entries such as name[n] into groups of the same n
707        dependent = {} # type: Dict[str, List[Parameter]]
708        dependent.update((p.id, []) for p in control)
709        for p in self.kernel_parameters:
710            if p.length_control is not None:
711                dependent[p.length_control].append(p)
712
713        # Gather entries such as name[4] into groups of the same length
714        fixed_length = {}  # type: Dict[int, List[Parameter]]
715        for p in self.kernel_parameters:
716            if p.length > 1 and p.length_control is None:
717                fixed_length.setdefault(p.length, []).append(p)
718
719        # Using the call_parameters table, we already have expanded forms
720        # for each of the vector parameters; put them in a lookup table
721        # Note: p.id and p.name are currently identical for the call parameters
722        expanded_pars = dict((p.id, p) for p in self.call_parameters)
723
724        def append_group(name):
725            """add the named parameter, and related magnetic parameters if any"""
726            result.append(expanded_pars[name])
727            if is2d:
728                for tag in '_M0', '_mtheta', '_mphi':
729                    if name+tag in expanded_pars:
730                        result.append(expanded_pars[name+tag])
731
732        # Gather the user parameters in order
733        result = control + self.common_parameters
734        for p in self.kernel_parameters:
735            if not is2d and p.type in ('orientation', 'magnetic'):
736                pass
737            elif p.is_control:
738                pass # already added
739            elif p.length_control is not None:
740                table = dependent.get(p.length_control, [])
741                if table:
742                    # look up length from incoming parameters
743                    table_length = int(pars.get(p.length_control, p.length))
744                    del dependent[p.length_control] # first entry seen
745                    for k in range(1, table_length+1):
746                        for entry in table:
747                            append_group(entry.id+str(k))
748                else:
749                    pass # already processed all entries
750            elif p.length > 1:
751                table = fixed_length.get(p.length, [])
752                if table:
753                    table_length = p.length
754                    del fixed_length[p.length]
755                    for k in range(1, table_length+1):
756                        for entry in table:
757                            append_group(entry.id+str(k))
758                else:
759                    pass # already processed all entries
760            else:
761                append_group(p.id)
762
763        if is2d and 'up_angle' in expanded_pars:
764            result.extend([
765                expanded_pars['up_frac_i'],
766                expanded_pars['up_frac_f'],
767                expanded_pars['up_angle'],
768            ])
769
770        return result
771
772def isstr(x):
773    # type: (Any) -> bool
774    """
775    Return True if the object is a string.
776    """
777    # TODO: 2-3 compatible tests for str, including unicode strings
778    return isinstance(x, str)
779
780
781#: Set of variables defined in the model that might contain C code
782C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code']
783
784def _find_source_lines(model_info, kernel_module):
785    # type: (ModelInfo, ModuleType) -> None
786    """
787    Identify the location of the C source inside the model definition file.
788
789    This code runs through the source of the kernel module looking for lines
790    that contain C code (because it is a c function definition).  Clearly
791    there are all sorts of reasons why this might not work (e.g., code
792    commented out in a triple-quoted line block, code built using string
793    concatenation, code defined in the branch of an 'if' block, code imported
794    from another file), but it should work properly in the 95% case, and for
795    the remainder, getting the incorrect line number will merely be
796    inconvenient.
797    """
798    # Only need line numbers if we are creating a C module and the C symbols
799    # are defined.
800    if (callable(model_info.Iq)
801            or not any(hasattr(model_info, s) for s in C_SYMBOLS)):
802        return
803
804    # load the module source if we can
805    try:
806        source = inspect.getsource(kernel_module)
807    except IOError:
808        return
809
810    # look for symbol at the start of the line
811    for lineno, line in enumerate(source.split('\n')):
812        for name in C_SYMBOLS:
813            if line.startswith(name):
814                # Add 1 since some compilers complain about "#line 0"
815                model_info.lineno[name] = lineno + 1
816                break
817
818def make_model_info(kernel_module):
819    # type: (module) -> ModelInfo
820    """
821    Extract the model definition from the loaded kernel module.
822
823    Fill in default values for parts of the module that are not provided.
824
825    Note: vectorized Iq and Iqac/Iqabc functions will be created for python
826    models when the model is first called, not when the model is loaded.
827    """
828    if hasattr(kernel_module, "model_info"):
829        # Custom sum/multi models
830        return kernel_module.model_info
831
832    info = ModelInfo()
833
834    # Build the parameter table
835    #print("make parameter table", kernel_module.parameters)
836    parameters = make_parameter_table(getattr(kernel_module, 'parameters', []))
837
838    # background defaults to zero for structure factor models
839    structure_factor = getattr(kernel_module, 'structure_factor', False)
840    if structure_factor:
841        parameters.set_zero_background()
842
843    # TODO: remove demo parameters
844    # The plots in the docs are generated from the model default values.
845    # Sascomp set parameters from the command line, and so doesn't need
846    # demo values for testing.
847    demo = expand_pars(parameters, getattr(kernel_module, 'demo', None))
848
849    filename = abspath(kernel_module.__file__).replace('.pyc', '.py')
850    kernel_id = splitext(basename(filename))[0]
851    name = getattr(kernel_module, 'name', None)
852    if name is None:
853        name = " ".join(w.capitalize() for w in kernel_id.split('_'))
854
855    info.id = kernel_id  # string used to load the kernel
856    info.filename = filename
857    info.name = name
858    info.title = getattr(kernel_module, 'title', name+" model")
859    info.description = getattr(kernel_module, 'description', 'no description')
860    info.parameters = parameters
861    info.demo = demo
862    info.composition = None
863    info.docs = kernel_module.__doc__
864    info.category = getattr(kernel_module, 'category', None)
865    info.structure_factor = getattr(kernel_module, 'structure_factor', False)
866    # TODO: find Fq by inspection
867    info.radius_effective_modes = getattr(kernel_module, 'radius_effective_modes', None)
868    info.have_Fq = getattr(kernel_module, 'have_Fq', False)
869    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
870    # Note: custom.load_custom_kernel_module assumes the C sources are defined
871    # by this attribute.
872    info.source = getattr(kernel_module, 'source', [])
873    info.c_code = getattr(kernel_module, 'c_code', None)
874    info.radius_effective = getattr(kernel_module, 'radius_effective', None)
875    # TODO: check the structure of the tests
876    info.tests = getattr(kernel_module, 'tests', [])
877    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore
878    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore
879    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore
880    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore
881    info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore
882    info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore
883    info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore
884    info.profile = getattr(kernel_module, 'profile', None) # type: ignore
885    info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore
886    # Default single and opencl to True for C models.  Python models have callable Iq.
887    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
888    info.single = getattr(kernel_module, 'single', not callable(info.Iq))
889    info.random = getattr(kernel_module, 'random', None)
890    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore
891
892    # Set control flag for explicitly set parameters, e.g., in the RPA model.
893    control = getattr(kernel_module, 'control', None)
894    if control is not None:
895        parameters[control].is_control = True
896
897    if callable(info.Iq) and parameters.has_2d:
898        raise ValueError("oriented python models not supported")
899
900    info.lineno = {}
901    _find_source_lines(info, kernel_module)
902    return info
903
904class ModelInfo(object):
905    """
906    Interpret the model definition file, categorizing the parameters.
907
908    The module can be loaded with a normal python import statement if you
909    know which module you need, or with __import__('sasmodels.model.'+name)
910    if the name is in a string.
911
912    The structure should be mostly static, other than the delayed definition
913    of *Iq*, *Iqac* and *Iqabc* if they need to be defined.
914    """
915    #: Full path to the file defining the kernel, if any.
916    filename = None         # type: Optional[str]
917    #: Id of the kernel used to load it from the filesystem.
918    id = None               # type: str
919    #: Display name of the model, which defaults to the model id but with
920    #: capitalization of the parts so for example core_shell defaults to
921    #: "Core Shell".
922    name = None             # type: str
923    #: Short description of the model.
924    title = None            # type: str
925    #: Long description of the model.
926    description = None      # type: str
927    #: Model parameter table. Parameters are defined using a list of parameter
928    #: definitions, each of which is contains parameter name, units,
929    #: default value, limits, type and description.  See :class:`Parameter`
930    #: for details on the individual parameters.  The parameters are gathered
931    #: into a :class:`ParameterTable`, which provides various views into the
932    #: parameter list.
933    parameters = None       # type: ParameterTable
934    #: Demo parameters as a *parameter:value* map used as the default values
935    #: for :mod:`compare`.  Any parameters not set in *demo* will use the
936    #: defaults from the parameter table.  That means no polydispersity, and
937    #: in the case of multiplicity models, a minimal model with no interesting
938    #: scattering.
939    demo = None             # type: Dict[str, float]
940    #: Composition is None if this is an independent model, or it is a
941    #: tuple with comoposition type ('product' or 'misture') and a list of
942    #: :class:`ModelInfo` blocks for the composed objects.  This allows us
943    #: to rebuild a complete mixture or product model from the info block.
944    #: *composition* is not given in the model definition file, but instead
945    #: arises when the model is constructed using names such as
946    #: *sphere*hardsphere* or *cylinder+sphere*.
947    composition = None      # type: Optional[Tuple[str, List[ModelInfo]]]
948    #: Different variants require different parameters.  In order to show
949    #: just the parameters needed for the variant selected by :attr:`control`,
950    #: you should provide a function *hidden(control) -> set(['a', 'b', ...])*
951    #: indicating which parameters need to be hidden.  For multiplicity
952    #: models, you need to use the complete name of the parameter, including
953    #: its number.  So for example, if variant "a" uses only *sld1* and *sld2*,
954    #: then *sld3*, *sld4* and *sld5* of multiplicity parameter *sld[5]*
955    #: should be in the hidden set.
956    hidden = None           # type: Optional[Callable[[int], Set[str]]]
957    #: Doc string from the top of the model file.  This should be formatted
958    #: using ReStructuredText format, with latex markup in ".. math"
959    #: environments, or in dollar signs.  This will be automatically
960    #: extracted to a .rst file by :func:`generate.make_docs`, then
961    #: converted to HTML or PDF by Sphinx.
962    docs = None             # type: str
963    #: Location of the model description in the documentation.  This takes the
964    #: form of "section" or "section:subsection".  So for example,
965    #: :ref:`porod` uses *category="shape-independent"* so it is in the
966    #: :ref:`shape-independent` section whereas
967    #: :ref:`capped-cylinder` uses: *category="shape:cylinder"*, which puts
968    #: it in the :ref:`shape-cylinder` section.
969    category = None         # type: Optional[str]
970    #: True if the model can be computed accurately with single precision.
971    #: This is True by default, but models such as :ref:`bcc-paracrystal` set
972    #: it to False because they require double precision calculations.
973    single = None           # type: bool
974    #: True if the model can be run as an opencl model.  If for some reason
975    #: the model cannot be run in opencl (e.g., because the model passes
976    #: functions by reference), then set this to false.
977    opencl = None           # type: bool
978    #: True if the model is a structure factor used to model the interaction
979    #: between form factor models.  This will default to False if it is not
980    #: provided in the file.
981    structure_factor = None # type: bool
982    #: True if the model defines an Fq function with signature
983    #: void Fq(double q, double *F1, double *F2, ...)
984    have_Fq = False
985    #: List of options for computing the effective radius of the shape,
986    #: or None if the model is not usable as a form factor model.
987    radius_effective_modes = None   # type: List[str]
988    #: List of C source files used to define the model.  The source files
989    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the
990    #: model defines orientation parameters. Files containing the most basic
991    #: functions must appear first in the list, followed by the files that
992    #: use those functions.
993    source = None           # type: List[str]
994    #: inline source code, added after all elements of source
995    c_code = None           # type: Optional[str]
996    #: Returns the form volume for python-based models.  Form volume is needed
997    #: for volume normalization in the polydispersity integral.  If no
998    #: parameters are *volume* parameters, then form volume is not needed.
999    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq`
1000    #: defined using a string containing C code), form_volume must also be
1001    #: C code, either defined as a string, or in the sources.
1002    form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]]
1003    #: Returns the shell volume for python-based models.  Form volume and
1004    #: shell volume are needed for volume normalization in the polydispersity
1005    #: integral and structure interactions for hollow shapes.  If no
1006    #: parameters are *volume* parameters, then shell volume is not needed.
1007    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq`
1008    #: defined using a string containing C code), shell_volume must also be
1009    #: C code, either defined as a string, or in the sources.
1010    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]]
1011    #: Computes the effective radius of the shape given the volume parameters.
1012    #: Only needed for models defined in python that can be used for
1013    #: monodisperse approximation for non-dilute solutions, P@S.  The first
1014    #: argument is the integer effective radius mode, with default 0.
1015    radius_effective = None  # type: Union[None, Callable[[int, np.ndarray], float]]
1016    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined
1017    #: by the parameter table.  *Iq* can be defined as a python function, or
1018    #: as a C function.  If it is defined in C, then set *Iq* to the body of
1019    #: the C function, including the return statement.  This function takes
1020    #: values for *q* and each of the parameters as separate *double* values
1021    #: (which may be converted to float or long double by sasmodels).  All
1022    #: source code files listed in :attr:`sources` will be loaded before the
1023    #: *Iq* function is defined.  If *Iq* is not present, then sources should
1024    #: define *static double Iq(double q, double a, double b, ...)* which
1025    #: will return *I(q, a, b, ...)*.  Multiplicity parameters are sent as
1026    #: pointers to doubles.  Constants in floating point expressions should
1027    #: include the decimal point. See :mod:`generate` for more details. If
1028    #: *have_Fq* is True, then Iq should return an interleaved array of
1029    #: $[\sum F(q_1), \sum F^2(q_1), \ldots, \sum F(q_n), \sum F^2(q_n)]$.
1030    Iq = None               # type: Union[None, str, Callable[[...], np.ndarray]]
1031    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
1032    Iqxy = None             # type: Union[None, str, Callable[[...], np.ndarray]]
1033    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
1034    Iqac = None             # type: Union[None, str, Callable[[...], np.ndarray]]
1035    #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`.
1036    Iqabc = None            # type: Union[None, str, Callable[[...], np.ndarray]]
1037    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
1038    Imagnetic = None        # type: Union[None, str, Callable[[...], np.ndarray]]
1039    #: Returns a model profile curve *x, y*.  If *profile* is defined, this
1040    #: curve will appear in response to the *Show* button in SasView.  Use
1041    #: :attr:`profile_axes` to set the axis labels.  Note that *y* values
1042    #: will be scaled by 1e6 before plotting.
1043    profile = None          # type: Optional[Callable[[np.ndarray], None]]
1044    #: Axis labels for the :attr:`profile` plot.  The default is *['x', 'y']*.
1045    #: Only the *x* component is used for now.
1046    profile_axes = None     # type: Tuple[str, str]
1047    #: Returns *sesans(z, a, b, ...)* for models which can directly compute
1048    #: the SESANS correlation function.  Note: not currently implemented.
1049    sesans = None           # type: Optional[Callable[[np.ndarray], np.ndarray]]
1050    #: Returns a random parameter set for the model
1051    random = None           # type: Optional[Callable[[], Dict[str, float]]]
1052    #: Line numbers for symbols defining C code
1053    lineno = None           # type: Dict[str, int]
1054    #: The set of tests that must pass.  The format of the tests is described
1055    #: in :mod:`model_test`.
1056    tests = None            # type: List[TestCondition]
1057
1058    def __init__(self):
1059        # type: () -> None
1060        pass
1061
1062    def get_hidden_parameters(self, control):
1063        """
1064        Returns the set of hidden parameters for the model.  *control* is the
1065        value of the control parameter.  Note that multiplicity models have
1066        an implicit control parameter, which is the parameter that controls
1067        the multiplicity.
1068        """
1069        if self.hidden is not None:
1070            hidden = self.hidden(control)
1071        else:
1072            controls = [p for p in self.parameters.kernel_parameters
1073                        if p.is_control]
1074            if len(controls) != 1:
1075                raise ValueError("more than one control parameter")
1076            hidden = set(p.id+str(k)
1077                         for p in self.parameters.kernel_parameters
1078                         for k in range(control+1, p.length+1)
1079                         if p.length > 1)
1080            for p in self.parameters.kernel_parameters:
1081                if p.length > 1 and p.type == "sld":
1082                    for k in range(control+1, p.length+1):
1083                        base = p.id+str(k)
1084                        hidden.update((base+"_M0", base+"_mtheta", base+"_mphi"))
1085        return hidden
Note: See TracBrowser for help on using the repository browser.