source: sasmodels/sasmodels/direct_model.py @ 65ceb7d

costrafo411
Last change on this file since 65ceb7d was 65ceb7d, checked in by Paul Kienzle <pkienzle@…>, 4 years ago

Merge branch 'master' into costrafo411

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