source: sasmodels/sasmodels/modelinfo.py @ d19962c

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

working vector parameter example using dll engine
  • Property mode set to 100644
File size: 18.4 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 = [1, len(choices)]
48    else:
49        choices = None
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
117class Parameter(object):
118    """
119    The available kernel parameters are defined as a list, with each parameter
120    defined as a sublist with the following elements:
121
122    *name* is the name that will be used in the call to the kernel
123    function and the name that will be displayed to the user.  Names
124    should be lower case, with words separated by underscore.  If
125    acronyms are used, the whole acronym should be upper case.
126
127    *units* should be one of *degrees* for angles, *Ang* for lengths,
128    *1e-6/Ang^2* for SLDs.
129
130    *default value* will be the initial value for  the model when it
131    is selected, or when an initial value is not otherwise specified.
132
133    *limits = [lb, ub]* are the hard limits on the parameter value, used to
134    limit the polydispersity density function.  In the fit, the parameter limits
135    given to the fit are the limits  on the central value of the parameter.
136    If there is polydispersity, it will evaluate parameter values outside
137    the fit limits, but not outside the hard limits specified in the model.
138    If there are no limits, use +/-inf imported from numpy.
139
140    *type* indicates how the parameter will be used.  "volume" parameters
141    will be used in all functions.  "orientation" parameters will be used
142    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in
143    *Imagnetic* only.  If *type* is the empty string, the parameter will
144    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters
145    can automatically be promoted to magnetic parameters, each of which
146    will have a magnitude and a direction, which may be different from
147    other sld parameters.
148
149    *description* is a short description of the parameter.  This will
150    be displayed in the parameter table and used as a tool tip for the
151    parameter value in the user interface.
152
153    Additional values can be set after the parameter is created:
154
155    *length* is the length of the field if it is a vector field
156    *length_control* is the parameter which sets the vector length
157    *is_control* is True if the parameter is a control parameter for a vector
158    *polydisperse* is true if the parameter accepts a polydispersity
159    *relative_pd* is true if that polydispersity is relative
160
161    In the usual process these values are set by :func:`make_parameter_table`
162    and :func:`parse_parameter` therein.
163    """
164    def __init__(self, name, units='', default=None, limits=(-np.inf, np.inf),
165                 type='', description=''):
166        self.id = name.split('[')[0].strip()
167        self.name = name
168        self.units = units
169        self.default = default
170        self.limits = limits
171        self.type = type
172        self.description = description
173        self.choices = None
174
175        # Length and length_control will be filled in once the complete
176        # parameter table is available.
177        self.length = 1
178        self.length_control = None
179        self.is_control = False
180
181        # TODO: need better control over whether a parameter is polydisperse
182        self.polydisperse = False
183        self.relative_pd = False
184
185    def as_definition(self):
186        """
187        Declare space for the variable in a parameter structure.
188
189        For example, the parameter thickness with length 3 will
190        return "double thickness[3];", with no spaces before and
191        no new line character afterward.
192        """
193        if self.length == 1:
194            return "double %s;"%self.id
195        else:
196            return "double %s[%d];"%(self.id, self.length)
197
198    def as_function_argument(self):
199        """
200        Declare the variable as a function argument.
201
202        For example, the parameter thickness with length 3 will
203        return "double *thickness", with no spaces before and
204        no comma afterward.
205        """
206        if self.length == 1:
207            return "double %s"%self.id
208        else:
209            return "double *%s"%self.id
210
211    def as_call_reference(self, prefix=""):
212        # Note: if the parameter is a struct type, then we will need to use
213        # &prefix+id.  For scalars and vectors we can just use prefix+id.
214        return prefix + self.id
215
216    def __str__(self):
217        return "<%s>"%self.name
218
219    def __repr__(self):
220        return "P<%s>"%self.name
221
222
223class ParameterTable(object):
224    """
225    ParameterTable manages the list of available parameters.
226
227    There are a couple of complications which mean that the list of parameters
228    for the kernel differs from the list of parameters that the user sees.
229
230    (1) Common parameters.  Scale and background are implicit to every model,
231    but are not passed to the kernel.
232
233    (2) Vector parameters.  Vector parameters are passed to the kernel as a
234    pointer to an array, e.g., thick[], but they are seen by the user as n
235    separate parameters thick1, thick2, ...
236
237    Therefore, the parameter table is organized by how it is expected to be
238    used. The following information is needed to set up the kernel functions:
239
240    * *kernel_parameters* is the list of parameters in the kernel parameter
241    table, with vector parameter p declared as p[].
242
243    * *iq_parameters* is the list of parameters to the Iq(q, ...) function,
244    with vector parameter p sent as p[].
245
246    * *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)
247    function, with vector parameter p sent as p[].
248
249    * *form_volume_parameters* is the list of parameters to the form_volume(...)
250    function, with vector parameter p sent as p[].
251
252    Problem details, which sets up the polydispersity loops, requires the
253    following:
254
255    * *theta_offset* is the offset of the theta parameter in the kernel parameter
256    table, with vector parameters counted as n individual parameters
257    p1, p2, ..., or offset is -1 if there is no theta parameter.
258
259    * *max_pd* is the maximum number of polydisperse parameters, with vector
260    parameters counted as n individual parameters p1, p2, ...  Note that
261    this number is limited to sasmodels.modelinfo.MAX_PD.
262
263    * *npars* is the total number of parameters to the kernel, with vector
264    parameters counted as n individual parameters p1, p2, ...
265
266    * *call_parameters* is the complete list of parameters to the kernel,
267    including scale and background, with vector parameters recorded as
268    individual parameters p1, p2, ...
269
270    * *active_1d* is the set of names that may be polydisperse for 1d data
271
272    * *active_2d* is the set of names that may be polydisperse for 2d data
273
274    User parameters are the set of parameters visible to the user, including
275    the scale and background parameters that the kernel does not see.  User
276    parameters don't use vector notation, and instead use p1, p2, ...
277
278    * *control_parameters* is the
279
280    """
281    # scale and background are implicit parameters
282    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS]
283
284    def __init__(self, parameters):
285        self.kernel_parameters = parameters
286        self._set_vector_lengths()
287        self._make_call_parameter_list()
288        self._categorize_parameters()
289        self._set_defaults()
290        #self._name_table= dict((p.id, p) for p in parameters)
291
292    def _set_vector_lengths(self):
293        # Sort out the length of the vector parameters such as thickness[n]
294        for p in self.kernel_parameters:
295            if p.length_control:
296                for ref in self.kernel_parameters:
297                    if ref.id == p.length_control:
298                        break
299                else:
300                    raise ValueError("no reference variable %r for %s"
301                                     % (p.length_control, p.name))
302                ref.is_control = True
303                low, high = ref.limits
304                if int(low) != low or int(high) != high or low<0 or high>20:
305                    raise ValueError("expected limits on %s to be within [0, 20]"
306                                     % ref.name)
307                p.length = high
308
309    def _set_defaults(self):
310        # Construct default values, including vector defaults
311        defaults = {}
312        for p in self.call_parameters:
313            if p.length == 1:
314                defaults[p.id] = p.default
315            else:
316                for k in range(1, p.length+1):
317                    defaults["%s%d"%(p.id, k)] = p.default
318        self.defaults = defaults
319
320    def _make_call_parameter_list(self):
321        full_list = self.COMMON[:]
322        for p in self.kernel_parameters:
323            if p.length == 1:
324                full_list.append(p)
325            else:
326                for k in range(1, p.length+1):
327                    pk = Parameter(p.id+str(k), p.units, p.default,
328                                   p.limits, p.type, p.description)
329                    pk.polydisperse = p.polydisperse
330                    pk.relative_pd = p.relative_pd
331                    full_list.append(pk)
332        self.call_parameters = full_list
333
334    """ # Suppress these for now until we see how they are used
335    def __getitem__(self, k):
336        if isinstance(k, (int, slice)):
337            return self.parameters[k]
338        else:
339            return self._name_table[k]
340
341    def __contains__(self, key):
342        return key in self._name_table
343
344    def __iter__(self):
345        return iter(self.expanded_parameters)
346    """
347
348    def _categorize_parameters(self):
349        # Set the kernel parameters.  Assumes background and scale are the
350        # first two parameters in the parameter list, but these are not sent
351        # to the underlying kernel functions.
352        self.iq_parameters = [p for p in self.kernel_parameters
353                              if p.type not in ('orientation', 'magnetic')]
354        self.iqxy_parameters = [p for p in self.kernel_parameters
355                                if p.type != 'magnetic']
356        self.form_volume_parameters = [p for p in self.kernel_parameters
357                                       if p.type == 'volume']
358
359        # Theta offset
360        offset = 0
361        for p in self.kernel_parameters:
362            if p.name == 'theta':
363                self.theta_offset = offset
364                break
365            offset += p.length
366        else:
367            self.theta_offset = -1
368
369        # number of polydisperse parameters
370        num_pd = sum(p.length for p in self.kernel_parameters if p.polydisperse)
371        # Don't use more polydisperse parameters than are available in the model
372        # Note: we can do polydispersity on arbitrary parameters, so it is not
373        # clear that this is a good idea; it does however make the poly_details
374        # code easier to write, so we will leave it in for now.
375        self.max_pd = min(num_pd, MAX_PD)
376
377        self.npars = sum(p.length for p in self.kernel_parameters)
378
379        # true if has 2D parameters
380        self.has_2d = any(p.type in ('orientation', 'magnetic')
381                          for p in self.kernel_parameters)
382
383        self.pd_1d = set(p.name for p in self.call_parameters
384                if p.polydisperse and p.type not in ('orientation', 'magnetic'))
385        self.pd_2d = set(p.name for p in self.call_parameters
386                         if p.polydisperse and p.type != 'magnetic')
387
388    def user_parameters(self, pars, is2d):
389        """
390        Return the list of parameters for the given data type.
391
392        Vector parameters are expanded as in place.  If multiple parameters
393        share the same vector length, then the parameters will be interleaved
394        in the result.  The control parameters come first.  For example,
395        if the parameter table is ordered as::
396
397            sld_core
398            sld_shell[num_shells]
399            sld_solvent
400            thickness[num_shells]
401            num_shells
402
403        and *pars[num_shells]=2* then the returned list will be::
404
405            num_shells
406            scale
407            background
408            sld_core
409            sld_shell1
410            thickness1
411            sld_shell2
412            thickness2
413            sld_solvent
414
415        Note that shell/thickness pairs are grouped together in the result
416        even though they were not grouped in the incoming table.  The control
417        parameter is always returned first since the GUI will want to set it
418        early, and rerender the table when it is changed.
419        """
420        control = [p for p in self.kernel_parameters if p.is_control]
421
422        # Gather entries such as name[n] into groups of the same n
423        dependent = dict((p.id, []) for p in control)
424        for p in self.kernel_parameters:
425            if p.length_control is not None:
426                dependent[p.length_control].append(p)
427
428        # Gather entries such as name[4] into groups of the same length
429        fixed = {}
430        for p in self.kernel_parameters:
431            if p.length > 1 and p.length_control is None:
432                fixed.setdefault(p.length, []).append(p)
433
434        # Using the call_parameters table, we already have expanded forms
435        # for each of the vector parameters; put them in a lookup table
436        expanded_pars = dict((p.name, p) for p in self.call_parameters)
437
438        # Gather the user parameters in order
439        result = control + self.COMMON
440        for p in self.kernel_parameters:
441            if not is2d and p.type in ('orientation', 'magnetic'):
442                pass
443            elif p.is_control:
444                pass # already added
445            elif p.length_control is not None:
446                table = dependent.get(p.length_control, [])
447                if table:
448                    # look up length from incoming parameters
449                    table_length = int(pars[p.length_control])
450                    del dependent[p.length_control] # first entry seen
451                    for k in range(1, table_length+1):
452                        for entry in table:
453                            result.append(expanded_pars[entry.id+str(k)])
454                else:
455                    pass # already processed all entries
456            elif p.length > 1:
457                table = fixed.get(p.length, [])
458                if table:
459                    table_length = p.length
460                    del fixed[p.length]
461                    for k in range(1, table_length+1):
462                        for entry in table:
463                            result.append(expanded_pars[entry.id+str(k)])
464                else:
465                    pass # already processed all entries
466            else:
467                result.append(p)
468
469        return result
470
471
472
473def isstr(x):
474    # TODO: 2-3 compatible tests for str, including unicode strings
475    return isinstance(x, str)
476
Note: See TracBrowser for help on using the repository browser.