source: sasmodels/sasmodels/modelinfo.py @ ce896fd

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

improved handling of vector parameters; remove compile errors from onion.c
  • Property mode set to 100644
File size: 20.1 KB
Line 
1
2import numpy as np
3
4# TODO: turn ModelInfo into a proper class
5ModelInfo = dict
6
7MAX_PD = 4
8
9COMMON_PARAMETERS = [
10    ["scale", "", 1, [0, np.inf], "", "Source intensity"],
11    ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"],
12]
13assert (len(COMMON_PARAMETERS) == 2
14        and COMMON_PARAMETERS[0][0]=="scale"
15        and COMMON_PARAMETERS[1][0]=="background"), "don't change common parameters"
16# assumptions about common parameters exist throughout the code, such as:
17# (1) kernel functions Iq, Iqxy, form_volume, ... don't see them
18# (2) kernel drivers assume scale is par[0] and background is par[1]
19# (3) mixture models drop the background on components and replace the scale
20#     with a scale that varies from [-inf, inf]
21# (4) product models drop the background and reassign scale
22# and maybe other places.
23# Note that scale and background cannot be coordinated parameters whose value
24# depends on the some polydisperse parameter with the current implementation
25
26def make_parameter_table(pars):
27    processed = []
28    for p in pars:
29        if not isinstance(p, list) or len(p) != 6:
30            raise ValueError("Parameter should be [name, units, default, limits, type, desc], but got %r"
31                             %str(p))
32        processed.append(parse_parameter(*p))
33    partable = ParameterTable(processed)
34    return partable
35
36def parse_parameter(name, units='', default=None,
37                    limits=(-np.inf, np.inf), type='', description=''):
38    # Parameter is a user facing class.  Do robust type checking.
39    if not isstr(name):
40        raise ValueError("expected string for parameter name %r"%name)
41    if not isstr(units):
42        raise ValueError("expected units to be a string for %s"%name)
43    # if limits is a list of strings, then this is a choice list
44    # field, and limits are 1 to length of string list
45    if isinstance(limits, list) and all(isstr(k) for k in limits):
46        choices = limits
47        limits = [0, len(choices)-1]
48    else:
49        choices = []
50    # TODO: maybe allow limits of None for (-inf, inf)
51    try:
52        low, high = limits
53        if not isinstance(low, (int, float)):
54            raise TypeError("low is not numeric")
55        if not isinstance(high, (int, float)):
56            raise TypeError("high is not numeric")
57        if low >= high:
58            raise ValueError("require low < high")
59    except:
60        raise ValueError("invalid limits %s for %s"%(limits, name))
61
62    if not isinstance(default, (int, float)):
63        raise ValueError("expected default %r to be a number for %s"
64                         % (default, name))
65    if default < low or default > high:
66        raise ValueError("default value %r not in range for %s"
67                         % (default, name))
68
69    if type not in ("volume", "orientation", "sld", "magnetic", ""):
70        raise ValueError("unexpected type %r for %s" % (type, name))
71
72    if not isstr(description):
73        raise ValueError("expected description to be a string")
74
75
76    # Parameter id for name[n] does not include [n]
77    if "[" in name:
78        if not name.endswith(']'):
79            raise ValueError("Expected name[len] for vector parameter %s"%name)
80        pid, ref = name[:-1].split('[', 1)
81        ref = ref.strip()
82    else:
83        pid, ref = name, None
84
85
86    # automatically identify sld types
87    if type=='' and (pid.startswith('sld') or pid.endswith('sld')):
88        type = 'sld'
89
90    # Check if using a vector definition, name[k], as the parameter name
91    if ref:
92        if ref == '':
93            raise ValueError("Need to specify vector length for %s"%name)
94        try:
95            length = int(ref)
96            control = None
97        except:
98            length = None
99            control = ref
100    else:
101        length = 1
102        control = None
103
104    # Build the parameter
105    parameter = Parameter(name=name, units=units, default=default,
106                          limits=limits, type=type, description=description)
107
108    # TODO: need better control over whether a parameter is polydisperse
109    parameter.polydisperse = type in ('orientation', 'volume')
110    parameter.relative_pd = type in ('volume')
111    parameter.choices = choices
112    parameter.length = length
113    parameter.length_control = control
114
115    return parameter
116
117
118def set_demo(model_info, demo):
119    """
120    Assign demo parameters to model_info['demo']
121
122    If demo is not defined, then use defaults.
123
124    If demo is defined, override defaults with value from demo.
125
126    If demo references vector fields, such as thickness[n], then support
127    different ways of assigning the demo values, including assigning a
128    specific value (e.g., thickness3=50.0), assigning a new value to all
129    (e.g., thickness=50.0) or assigning values using list notation.
130    """
131    partable = model_info['parameters']
132    if demo is None:
133        result = partable.defaults
134    else:
135        pars = dict((p.id, p) for p in partable.kernel_parameters)
136        result = partable.defaults.copy()
137        vectors = dict((name,value) for name,value in demo.items()
138                       if name in pars and pars[name].length > 1)
139        if vectors:
140            scalars = dict((name, value) for name, value in demo.items()
141                           if name not in pars or pars[name].length == 1)
142            for name, value in vectors.items():
143                if np.isscalar(value):
144                    # support for the form
145                    #    demo(thickness=0, thickness2=50)
146                    for k in pars[name].length:
147                        key = name+str(k)
148                        if key not in scalars:
149                            scalars[key] = vectors
150                else:
151                    # supoprt for the form
152                    #    demo(thickness=[20,10,3])
153                    for (k,v) in enumerate(value):
154                        scalars[name+str(k)] = v
155            result.update(scalars)
156        else:
157            result.update(demo)
158
159    model_info['demo'] = result
160
161class Parameter(object):
162    """
163    The available kernel parameters are defined as a list, with each parameter
164    defined as a sublist with the following elements:
165
166    *name* is the name that will be used in the call to the kernel
167    function and the name that will be displayed to the user.  Names
168    should be lower case, with words separated by underscore.  If
169    acronyms are used, the whole acronym should be upper case.
170
171    *units* should be one of *degrees* for angles, *Ang* for lengths,
172    *1e-6/Ang^2* for SLDs.
173
174    *default value* will be the initial value for  the model when it
175    is selected, or when an initial value is not otherwise specified.
176
177    *limits = [lb, ub]* are the hard limits on the parameter value, used to
178    limit the polydispersity density function.  In the fit, the parameter limits
179    given to the fit are the limits  on the central value of the parameter.
180    If there is polydispersity, it will evaluate parameter values outside
181    the fit limits, but not outside the hard limits specified in the model.
182    If there are no limits, use +/-inf imported from numpy.
183
184    *type* indicates how the parameter will be used.  "volume" parameters
185    will be used in all functions.  "orientation" parameters will be used
186    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
187    *Imagnetic* only.  If *type* is the empty string, the parameter will
188    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters
189    can automatically be promoted to magnetic parameters, each of which
190    will have a magnitude and a direction, which may be different from
191    other sld parameters.
192
193    *description* is a short description of the parameter.  This will
194    be displayed in the parameter table and used as a tool tip for the
195    parameter value in the user interface.
196
197    Additional values can be set after the parameter is created:
198
199    *length* is the length of the field if it is a vector field
200    *length_control* is the parameter which sets the vector length
201    *is_control* is True if the parameter is a control parameter for a vector
202    *polydisperse* is true if the parameter accepts a polydispersity
203    *relative_pd* is true if that polydispersity is relative
204
205    In the usual process these values are set by :func:`make_parameter_table`
206    and :func:`parse_parameter` therein.
207    """
208    def __init__(self, name, units='', default=None, limits=(-np.inf, np.inf),
209                 type='', description=''):
210        self.id = name.split('[')[0].strip()
211        self.name = name
212        self.units = units
213        self.default = default
214        self.limits = limits
215        self.type = type
216        self.description = description
217
218        # Length and length_control will be filled in once the complete
219        # parameter table is available.
220        self.length = 1
221        self.length_control = None
222        self.is_control = False
223
224        # TODO: need better control over whether a parameter is polydisperse
225        self.polydisperse = False
226        self.relative_pd = False
227
228        # choices are also set externally.
229        self.choices = []
230
231    def as_definition(self):
232        """
233        Declare space for the variable in a parameter structure.
234
235        For example, the parameter thickness with length 3 will
236        return "double thickness[3];", with no spaces before and
237        no new line character afterward.
238        """
239        if self.length == 1:
240            return "double %s;"%self.id
241        else:
242            return "double %s[%d];"%(self.id, self.length)
243
244    def as_function_argument(self):
245        """
246        Declare the variable as a function argument.
247
248        For example, the parameter thickness with length 3 will
249        return "double *thickness", with no spaces before and
250        no comma afterward.
251        """
252        if self.length == 1:
253            return "double %s"%self.id
254        else:
255            return "double *%s"%self.id
256
257    def as_call_reference(self, prefix=""):
258        # Note: if the parameter is a struct type, then we will need to use
259        # &prefix+id.  For scalars and vectors we can just use prefix+id.
260        return prefix + self.id
261
262    def __str__(self):
263        return "<%s>"%self.name
264
265    def __repr__(self):
266        return "P<%s>"%self.name
267
268
269class ParameterTable(object):
270    """
271    ParameterTable manages the list of available parameters.
272
273    There are a couple of complications which mean that the list of parameters
274    for the kernel differs from the list of parameters that the user sees.
275
276    (1) Common parameters.  Scale and background are implicit to every model,
277    but are not passed to the kernel.
278
279    (2) Vector parameters.  Vector parameters are passed to the kernel as a
280    pointer to an array, e.g., thick[], but they are seen by the user as n
281    separate parameters thick1, thick2, ...
282
283    Therefore, the parameter table is organized by how it is expected to be
284    used. The following information is needed to set up the kernel functions:
285
286    * *kernel_parameters* is the list of parameters in the kernel parameter
287    table, with vector parameter p declared as p[].
288
289    * *iq_parameters* is the list of parameters to the Iq(q, ...) function,
290    with vector parameter p sent as p[].
291
292    * *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)
293    function, with vector parameter p sent as p[].
294
295    * *form_volume_parameters* is the list of parameters to the form_volume(...)
296    function, with vector parameter p sent as p[].
297
298    Problem details, which sets up the polydispersity loops, requires the
299    following:
300
301    * *theta_offset* is the offset of the theta parameter in the kernel parameter
302    table, with vector parameters counted as n individual parameters
303    p1, p2, ..., or offset is -1 if there is no theta parameter.
304
305    * *max_pd* is the maximum number of polydisperse parameters, with vector
306    parameters counted as n individual parameters p1, p2, ...  Note that
307    this number is limited to sasmodels.modelinfo.MAX_PD.
308
309    * *npars* is the total number of parameters to the kernel, with vector
310    parameters counted as n individual parameters p1, p2, ...
311
312    * *call_parameters* is the complete list of parameters to the kernel,
313    including scale and background, with vector parameters recorded as
314    individual parameters p1, p2, ...
315
316    * *active_1d* is the set of names that may be polydisperse for 1d data
317
318    * *active_2d* is the set of names that may be polydisperse for 2d data
319
320    User parameters are the set of parameters visible to the user, including
321    the scale and background parameters that the kernel does not see.  User
322    parameters don't use vector notation, and instead use p1, p2, ...
323
324    * *control_parameters* is the
325
326    """
327    # scale and background are implicit parameters
328    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS]
329
330    def __init__(self, parameters):
331        self.kernel_parameters = parameters
332        self._set_vector_lengths()
333        self._make_call_parameter_list()
334        self._categorize_parameters()
335        self._set_defaults()
336        #self._name_table= dict((p.id, p) for p in parameters)
337
338    def _set_vector_lengths(self):
339        # Sort out the length of the vector parameters such as thickness[n]
340        for p in self.kernel_parameters:
341            if p.length_control:
342                for ref in self.kernel_parameters:
343                    if ref.id == p.length_control:
344                        break
345                else:
346                    raise ValueError("no reference variable %r for %s"
347                                     % (p.length_control, p.name))
348                ref.is_control = True
349                low, high = ref.limits
350                if int(low) != low or int(high) != high or low<0 or high>20:
351                    raise ValueError("expected limits on %s to be within [0, 20]"
352                                     % ref.name)
353                p.length = high
354
355    def _set_defaults(self):
356        # Construct default values, including vector defaults
357        defaults = {}
358        for p in self.call_parameters:
359            if p.length == 1:
360                defaults[p.id] = p.default
361            else:
362                for k in range(1, p.length+1):
363                    defaults["%s%d"%(p.id, k)] = p.default
364        self.defaults = defaults
365
366    def _make_call_parameter_list(self):
367        full_list = self.COMMON[:]
368        for p in self.kernel_parameters:
369            if p.length == 1:
370                full_list.append(p)
371            else:
372                for k in range(1, p.length+1):
373                    pk = Parameter(p.id+str(k), p.units, p.default,
374                                   p.limits, p.type, p.description)
375                    pk.polydisperse = p.polydisperse
376                    pk.relative_pd = p.relative_pd
377                    full_list.append(pk)
378        self.call_parameters = full_list
379
380    """ # Suppress these for now until we see how they are used
381    def __getitem__(self, k):
382        if isinstance(k, (int, slice)):
383            return self.parameters[k]
384        else:
385            return self._name_table[k]
386
387    def __contains__(self, key):
388        return key in self._name_table
389
390    def __iter__(self):
391        return iter(self.expanded_parameters)
392    """
393
394    def _categorize_parameters(self):
395        # Set the kernel parameters.  Assumes background and scale are the
396        # first two parameters in the parameter list, but these are not sent
397        # to the underlying kernel functions.
398        self.iq_parameters = [p for p in self.kernel_parameters
399                              if p.type not in ('orientation', 'magnetic')]
400        self.iqxy_parameters = [p for p in self.kernel_parameters
401                                if p.type != 'magnetic']
402        self.form_volume_parameters = [p for p in self.kernel_parameters
403                                       if p.type == 'volume']
404
405        # Theta offset
406        offset = 0
407        for p in self.kernel_parameters:
408            if p.name == 'theta':
409                self.theta_offset = offset
410                break
411            offset += p.length
412        else:
413            self.theta_offset = -1
414
415        # number of polydisperse parameters
416        num_pd = sum(p.length for p in self.kernel_parameters if p.polydisperse)
417        # Don't use more polydisperse parameters than are available in the model
418        # Note: we can do polydispersity on arbitrary parameters, so it is not
419        # clear that this is a good idea; it does however make the poly_details
420        # code easier to write, so we will leave it in for now.
421        self.max_pd = min(num_pd, MAX_PD)
422
423        self.npars = sum(p.length for p in self.kernel_parameters)
424
425        # true if has 2D parameters
426        self.has_2d = any(p.type in ('orientation', 'magnetic')
427                          for p in self.kernel_parameters)
428
429        self.pd_1d = set(p.name for p in self.call_parameters
430                if p.polydisperse and p.type not in ('orientation', 'magnetic'))
431        self.pd_2d = set(p.name for p in self.call_parameters
432                         if p.polydisperse and p.type != 'magnetic')
433
434    def user_parameters(self, pars, is2d):
435        """
436        Return the list of parameters for the given data type.
437
438        Vector parameters are expanded as in place.  If multiple parameters
439        share the same vector length, then the parameters will be interleaved
440        in the result.  The control parameters come first.  For example,
441        if the parameter table is ordered as::
442
443            sld_core
444            sld_shell[num_shells]
445            sld_solvent
446            thickness[num_shells]
447            num_shells
448
449        and *pars[num_shells]=2* then the returned list will be::
450
451            num_shells
452            scale
453            background
454            sld_core
455            sld_shell1
456            thickness1
457            sld_shell2
458            thickness2
459            sld_solvent
460
461        Note that shell/thickness pairs are grouped together in the result
462        even though they were not grouped in the incoming table.  The control
463        parameter is always returned first since the GUI will want to set it
464        early, and rerender the table when it is changed.
465        """
466        control = [p for p in self.kernel_parameters if p.is_control]
467
468        # Gather entries such as name[n] into groups of the same n
469        dependent = dict((p.id, []) for p in control)
470        for p in self.kernel_parameters:
471            if p.length_control is not None:
472                dependent[p.length_control].append(p)
473
474        # Gather entries such as name[4] into groups of the same length
475        fixed = {}
476        for p in self.kernel_parameters:
477            if p.length > 1 and p.length_control is None:
478                fixed.setdefault(p.length, []).append(p)
479
480        # Using the call_parameters table, we already have expanded forms
481        # for each of the vector parameters; put them in a lookup table
482        expanded_pars = dict((p.name, p) for p in self.call_parameters)
483
484        # Gather the user parameters in order
485        result = control + self.COMMON
486        for p in self.kernel_parameters:
487            if not is2d and p.type in ('orientation', 'magnetic'):
488                pass
489            elif p.is_control:
490                pass # already added
491            elif p.length_control is not None:
492                table = dependent.get(p.length_control, [])
493                if table:
494                    # look up length from incoming parameters
495                    table_length = int(pars[p.length_control])
496                    del dependent[p.length_control] # first entry seen
497                    for k in range(1, table_length+1):
498                        for entry in table:
499                            result.append(expanded_pars[entry.id+str(k)])
500                else:
501                    pass # already processed all entries
502            elif p.length > 1:
503                table = fixed.get(p.length, [])
504                if table:
505                    table_length = p.length
506                    del fixed[p.length]
507                    for k in range(1, table_length+1):
508                        for entry in table:
509                            result.append(expanded_pars[entry.id+str(k)])
510                else:
511                    pass # already processed all entries
512            else:
513                result.append(p)
514
515        return result
516
517
518
519def isstr(x):
520    # TODO: 2-3 compatible tests for str, including unicode strings
521    return isinstance(x, str)
522
Note: See TracBrowser for help on using the repository browser.