source: sasmodels/sasmodels/direct_model.py @ a839b22

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

restore working sesans example using direct model

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