source: sasmodels/sasmodels/direct_model.py @ 3c24ccd

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3c24ccd was 3c24ccd, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

add -weights option to sascomp to show dispersity distribution

  • Property mode set to 100644
File size: 17.0 KB
Line 
1"""
2Class interface to the model calculator.
3
4Calling a model is somewhat non-trivial since the functions called depend
5on the data type.  For 1D data the *Iq* kernel needs to be called, for
62D data the *Iqxy* kernel needs to be called, and for SESANS data the
7*Iq* kernel needs to be called followed by a Hankel transform.  Before
8the kernel is called an appropriate *q* calculation vector needs to be
9constructed.  This is not the simple *q* vector where you have measured
10the data since the resolution calculation will require values beyond the
11range of the measured data.  After the calculation the resolution calculator
12must be called to return the predicted value for each measured data point.
13
14:class:`DirectModel` is a callable object that takes *parameter=value*
15keyword arguments and returns the appropriate theory values for the data.
16
17:class:`DataMixin` does the real work of interpreting the data and calling
18the model calculator.  This is used by :class:`DirectModel`, which uses
19direct parameter values and by :class:`bumps_model.Experiment` which wraps
20the parameter values in boxes so that the user can set fitting ranges, etc.
21on the individual parameters and send the model to the Bumps optimizers.
22"""
23from __future__ import print_function
24
25import numpy as np  # type: ignore
26
27# TODO: fix sesans module
28from . import sesans  # type: ignore
29from . import weights
30from . import resolution
31from . import resolution2d
32from .details import make_kernel_args, dispersion_mesh
33
34try:
35    from typing import Optional, Dict, Tuple
36except ImportError:
37    pass
38else:
39    from .data import Data
40    from .kernel import Kernel, KernelModel
41    from .modelinfo import Parameter, ParameterSet
42
43def call_kernel(calculator, pars, cutoff=0., mono=False):
44    # type: (Kernel, ParameterSet, float, bool) -> np.ndarray
45    """
46    Call *kernel* returned from *model.make_kernel* with parameters *pars*.
47
48    *cutoff* is the limiting value for the product of dispersion weights used
49    to perform the multidimensional dispersion calculation more quickly at a
50    slight cost to accuracy. The default value of *cutoff=0* integrates over
51    the entire dispersion cube.  Using *cutoff=1e-5* can be 50% faster, but
52    with an error of about 1%, which is usually less than the measurement
53    uncertainty.
54
55    *mono* is True if polydispersity should be set to none on all parameters.
56    """
57    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono)
58    call_details, values, is_magnetic = make_kernel_args(calculator, mesh)
59    #print("values:", values)
60    return calculator(call_details, values, cutoff, is_magnetic)
61
62def call_ER(model_info, pars):
63    # type: (ModelInfo, ParameterSet) -> float
64    """
65    Call the model ER function using *values*.
66
67    *model_info* is either *model.info* if you have a loaded model,
68    or *kernel.info* if you have a model kernel prepared for evaluation.
69    """
70    if model_info.ER is None:
71        return 1.0
72    elif not model_info.parameters.form_volume_parameters:
73        # handle the case where ER is provided but model is not polydisperse
74        return model_info.ER()
75    else:
76        value, weight = _vol_pars(model_info, pars)
77        individual_radii = model_info.ER(*value)
78        return np.sum(weight*individual_radii) / np.sum(weight)
79
80
81def call_VR(model_info, pars):
82    # type: (ModelInfo, ParameterSet) -> float
83    """
84    Call the model VR function using *pars*.
85
86    *model_info* is either *model.info* if you have a loaded model,
87    or *kernel.info* if you have a model kernel prepared for evaluation.
88    """
89    if model_info.VR is None:
90        return 1.0
91    elif not model_info.parameters.form_volume_parameters:
92        # handle the case where ER is provided but model is not polydisperse
93        return model_info.VR()
94    else:
95        value, weight = _vol_pars(model_info, pars)
96        whole, part = model_info.VR(*value)
97        return np.sum(weight*part)/np.sum(weight*whole)
98
99
100def call_profile(model_info, **pars):
101    # type: (ModelInfo, ...) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]]
102    """
103    Returns the profile *x, y, (xlabel, ylabel)* representing the model.
104    """
105    args = {}
106    for p in model_info.parameters.kernel_parameters:
107        if p.length > 1:
108            value = np.array([pars.get(p.id+str(j), p.default)
109                              for j in range(1, p.length+1)])
110        else:
111            value = pars.get(p.id, p.default)
112        args[p.id] = value
113    x, y = model_info.profile(**args)
114    return x, y, model_info.profile_axes
115
116def get_mesh(model_info, values, dim='1d', mono=False):
117    # type: (ModelInfo, Dict[str, float], str, bool) -> List[Tuple[float, np.ndarray, np.ndarry]]
118    """
119    Retrieve the dispersity mesh described by the parameter set.
120
121    Returns a list of *(value, dispersity, weights)* with one tuple for each
122    parameter in the model call parameters.  Inactive parameters return the
123    default value with a weight of 1.0.
124    """
125    parameters = model_info.parameters
126    if mono:
127        active = lambda name: False
128    elif dim == '1d':
129        active = lambda name: name in parameters.pd_1d
130    elif dim == '2d':
131        active = lambda name: name in parameters.pd_2d
132    else:
133        active = lambda name: True
134
135    #print("pars",[p.id for p in parameters.call_parameters])
136    mesh = [_get_par_weights(p, values, active(p.name))
137            for p in parameters.call_parameters]
138    return mesh
139
140
141def _get_par_weights(parameter, values, active=True):
142    # type: (Parameter, Dict[str, float]) -> Tuple[float, np.ndarray, np.ndarray]
143    """
144    Generate the distribution for parameter *name* given the parameter values
145    in *pars*.
146
147    Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma"
148    from the *pars* dictionary for parameter value and parameter dispersion.
149    """
150    value = float(values.get(parameter.name, parameter.default))
151    npts = values.get(parameter.name+'_pd_n', 0)
152    width = values.get(parameter.name+'_pd', 0.0)
153    relative = parameter.relative_pd
154    if npts == 0 or width == 0 or not active:
155        # Note: orientation parameters have the viewing angle as the parameter
156        # value and the jitter in the distribution, so be sure to set the
157        # empty pd for orientation parameters to 0.
158        pd = [value if relative else 0.0], [1.0]
159    else:
160        limits = parameter.limits
161        disperser = values.get(parameter.name+'_pd_type', 'gaussian')
162        nsigma = values.get(parameter.name+'_pd_nsigma', 3.0)
163        pd = weights.get_weights(disperser, npts, width, nsigma,
164                                value, limits, relative)
165    return value, pd[0], pd[1]
166
167
168def _vol_pars(model_info, pars):
169    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray]
170    vol_pars = [get_weights(p, pars)
171                for p in model_info.parameters.call_parameters
172                if p.type == 'volume']
173    #import pylab; pylab.plot(vol_pars[0][0],vol_pars[0][1]); pylab.show()
174    value, weight = dispersion_mesh(model_info, vol_pars)
175    return value, weight
176
177
178class DataMixin(object):
179    """
180    DataMixin captures the common aspects of evaluating a SAS model for a
181    particular data set, including calculating Iq and evaluating the
182    resolution function.  It is used in particular by :class:`DirectModel`,
183    which evaluates a SAS model parameters as key word arguments to the
184    calculator method, and by :class:`bumps_model.Experiment`, which wraps the
185    model and data for use with the Bumps fitting engine.  It is not
186    currently used by :class:`sasview_model.SasviewModel` since this will
187    require a number of changes to SasView before we can do it.
188
189    :meth:`_interpret_data` initializes the data structures necessary
190    to manage the calculations.  This sets attributes in the child class
191    such as *data_type* and *resolution*.
192
193    :meth:`_calc_theory` evaluates the model at the given control values.
194
195    :meth:`_set_data` sets the intensity data in the data object,
196    possibly with random noise added.  This is useful for simulating a
197    dataset with the results from :meth:`_calc_theory`.
198    """
199    def _interpret_data(self, data, model):
200        # type: (Data, KernelModel) -> None
201        # pylint: disable=attribute-defined-outside-init
202
203        self._data = data
204        self._model = model
205
206        # interpret data
207        if hasattr(data, 'isSesans') and data.isSesans:
208            self.data_type = 'sesans'
209        elif hasattr(data, 'qx_data'):
210            self.data_type = 'Iqxy'
211        elif getattr(data, 'oriented', False):
212            self.data_type = 'Iq-oriented'
213        else:
214            self.data_type = 'Iq'
215
216        if self.data_type == 'sesans':
217            q = sesans.make_q(data.sample.zacceptance, data.Rmax)
218            index = slice(None, None)
219            res = None
220            if data.y is not None:
221                Iq, dIq = data.y, data.dy
222            else:
223                Iq, dIq = None, None
224            #self._theory = np.zeros_like(q)
225            q_vectors = [q]
226            q_mono = sesans.make_all_q(data)
227        elif self.data_type == 'Iqxy':
228            #if not model.info.parameters.has_2d:
229            #    raise ValueError("not 2D without orientation or magnetic parameters")
230            q = np.sqrt(data.qx_data**2 + data.qy_data**2)
231            qmin = getattr(data, 'qmin', 1e-16)
232            qmax = getattr(data, 'qmax', np.inf)
233            accuracy = getattr(data, 'accuracy', 'Low')
234            index = ~data.mask & (q >= qmin) & (q <= qmax)
235            if data.data is not None:
236                index &= ~np.isnan(data.data)
237                Iq = data.data[index]
238                dIq = data.err_data[index]
239            else:
240                Iq, dIq = None, None
241            res = resolution2d.Pinhole2D(data=data, index=index,
242                                         nsigma=3.0, accuracy=accuracy)
243            #self._theory = np.zeros_like(self.Iq)
244            q_vectors = res.q_calc
245            q_mono = []
246        elif self.data_type == 'Iq':
247            index = (data.x >= data.qmin) & (data.x <= data.qmax)
248            if data.y is not None:
249                index &= ~np.isnan(data.y)
250                Iq = data.y[index]
251                dIq = data.dy[index]
252            else:
253                Iq, dIq = None, None
254            if getattr(data, 'dx', None) is not None:
255                q, dq = data.x[index], data.dx[index]
256                if (dq > 0).any():
257                    res = resolution.Pinhole1D(q, dq)
258                else:
259                    res = resolution.Perfect1D(q)
260            elif (getattr(data, 'dxl', None) is not None
261                  and getattr(data, 'dxw', None) is not None):
262                res = resolution.Slit1D(data.x[index],
263                                        qx_width=data.dxl[index],
264                                        qy_width=data.dxw[index])
265            else:
266                res = resolution.Perfect1D(data.x[index])
267
268            #self._theory = np.zeros_like(self.Iq)
269            q_vectors = [res.q_calc]
270            q_mono = []
271        elif self.data_type == 'Iq-oriented':
272            index = (data.x >= data.qmin) & (data.x <= data.qmax)
273            if data.y is not None:
274                index &= ~np.isnan(data.y)
275                Iq = data.y[index]
276                dIq = data.dy[index]
277            else:
278                Iq, dIq = None, None
279            if (getattr(data, 'dxl', None) is None
280                    or getattr(data, 'dxw', None) is None):
281                raise ValueError("oriented sample with 1D data needs slit resolution")
282
283            res = resolution2d.Slit2D(data.x[index],
284                                      qx_width=data.dxw[index],
285                                      qy_width=data.dxl[index])
286            q_vectors = res.q_calc
287            q_mono = []
288        else:
289            raise ValueError("Unknown data type") # never gets here
290
291        # Remember function inputs so we can delay loading the function and
292        # so we can save/restore state
293        self._kernel_inputs = q_vectors
294        self._kernel_mono_inputs = q_mono
295        self._kernel = None
296        self.Iq, self.dIq, self.index = Iq, dIq, index
297        self.resolution = res
298
299    def _set_data(self, Iq, noise=None):
300        # type: (np.ndarray, Optional[float]) -> None
301        # pylint: disable=attribute-defined-outside-init
302        if noise is not None:
303            self.dIq = Iq*noise*0.01
304        dy = self.dIq
305        y = Iq + np.random.randn(*dy.shape) * dy
306        self.Iq = y
307        if self.data_type in ('Iq', 'Iq-oriented'):
308            if self._data.y is None:
309                self._data.y = np.empty(len(self._data.x), 'd')
310            if self._data.dy is None:
311                self._data.dy = np.empty(len(self._data.x), 'd')
312            self._data.dy[self.index] = dy
313            self._data.y[self.index] = y
314        elif self.data_type == 'Iqxy':
315            if self._data.data is None:
316                self._data.data = np.empty_like(self._data.qx_data, 'd')
317            if self._data.err_data is None:
318                self._data.err_data = np.empty_like(self._data.qx_data, 'd')
319            self._data.data[self.index] = y
320            self._data.err_data[self.index] = dy
321        elif self.data_type == 'sesans':
322            if self._data.y is None:
323                self._data.y = np.empty(len(self._data.x), 'd')
324            self._data.y[self.index] = y
325        else:
326            raise ValueError("Unknown model")
327
328    def _calc_theory(self, pars, cutoff=0.0):
329        # type: (ParameterSet, float) -> np.ndarray
330        if self._kernel is None:
331            self._kernel = self._model.make_kernel(self._kernel_inputs)
332            self._kernel_mono = (
333                self._model.make_kernel(self._kernel_mono_inputs)
334                if self._kernel_mono_inputs else None)
335
336        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
337        # Storing the calculated Iq values so that they can be plotted.
338        # Only applies to oriented USANS data for now.
339        # TODO: extend plotting of calculate Iq to other measurement types
340        # TODO: refactor so we don't store the result in the model
341        self.Iq_calc = Iq_calc
342        if self.data_type == 'sesans':
343            Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True)
344                       if self._kernel_mono_inputs else None)
345            result = sesans.transform(self._data,
346                                      self._kernel_inputs[0], Iq_calc,
347                                      self._kernel_mono_inputs, Iq_mono)
348        else:
349            result = self.resolution.apply(Iq_calc)
350            if hasattr(self.resolution, 'nx'):
351                self.Iq_calc = (
352                    self.resolution.qx_calc, self.resolution.qy_calc,
353                    np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx))
354                )
355        return result
356
357
358class DirectModel(DataMixin):
359    """
360    Create a calculator object for a model.
361
362    *data* is 1D SAS, 2D SAS or SESANS data
363
364    *model* is a model calculator return from :func:`generate.load_model`
365
366    *cutoff* is the polydispersity weight cutoff.
367    """
368    def __init__(self, data, model, cutoff=1e-5):
369        # type: (Data, KernelModel, float) -> None
370        self.model = model
371        self.cutoff = cutoff
372        # Note: _interpret_data defines the model attributes
373        self._interpret_data(data, model)
374
375    def __call__(self, **pars):
376        # type: (**float) -> np.ndarray
377        return self._calc_theory(pars, cutoff=self.cutoff)
378
379    def simulate_data(self, noise=None, **pars):
380        # type: (Optional[float], **float) -> None
381        """
382        Generate simulated data for the model.
383        """
384        Iq = self.__call__(**pars)
385        self._set_data(Iq, noise=noise)
386
387    def profile(self, **pars):
388        # type: (**float) -> None
389        """
390        Generate a plottable profile.
391        """
392        return call_profile(self.model.info, **pars)
393
394def main():
395    # type: () -> None
396    """
397    Program to evaluate a particular model at a set of q values.
398    """
399    import sys
400    from .data import empty_data1D, empty_data2D
401    from .core import load_model_info, build_model
402
403    if len(sys.argv) < 3:
404        print("usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ...")
405        sys.exit(1)
406    model_name = sys.argv[1]
407    call = sys.argv[2].upper()
408    if call != "ER_VR":
409        try:
410            values = [float(v) for v in call.split(',')]
411        except Exception:
412            values = []
413        if len(values) == 1:
414            q, = values
415            data = empty_data1D([q])
416        elif len(values) == 2:
417            qx, qy = values
418            data = empty_data2D([qx], [qy])
419        else:
420            print("use q or qx,qy or ER or VR")
421            sys.exit(1)
422    else:
423        data = empty_data1D([0.001])  # Data not used in ER/VR
424
425    model_info = load_model_info(model_name)
426    model = build_model(model_info)
427    calculator = DirectModel(data, model)
428    pars = dict((k, (float(v) if not k.endswith("_pd_type") else v))
429                for pair in sys.argv[3:]
430                for k, v in [pair.split('=')])
431    if call == "ER_VR":
432        ER = call_ER(model_info, pars)
433        VR = call_VR(model_info, pars)
434        print(ER, VR)
435    else:
436        Iq = calculator(**pars)
437        print(Iq[0])
438
439if __name__ == "__main__":
440    main()
Note: See TracBrowser for help on using the repository browser.