source: sasmodels/sasmodels/direct_model.py @ 2cdc35b

costrafo411
Last change on this file since 2cdc35b was 2cdc35b, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

provide working oriented/unoriented sesans examples

  • Property mode set to 100644
File size: 15.6 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            self._data.dy[self.index] = dy
299            self._data.y[self.index] = y
300        elif self.data_type == 'Iqxy':
301            self._data.data[self.index] = y
302        elif self.data_type == 'sesans':
303            self._data.y[self.index] = y
304        else:
305            raise ValueError("Unknown model")
306
307    def _calc_theory(self, pars, cutoff=0.0):
308        # type: (ParameterSet, float) -> np.ndarray
309        if self._kernel is None:
310            self._kernel = self._model.make_kernel(self._kernel_inputs)
311
312        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
313        # Storing the calculated Iq values so that they can be plotted.
314        # Only applies to oriented USANS data for now.
315        # TODO: extend plotting of calculate Iq to other measurement types
316        # TODO: refactor so we don't store the result in the model
317        self.Iq_calc = None
318        result = self.resolution.apply(Iq_calc)
319        if hasattr(self.resolution, 'nx'):
320            self.Iq_calc = (
321                self.resolution.qx_calc, self.resolution.qy_calc,
322                np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx))
323            )
324        return result
325
326
327class DirectModel(DataMixin):
328    """
329    Create a calculator object for a model.
330
331    *data* is 1D SAS, 2D SAS or SESANS data
332
333    *model* is a model calculator return from :func:`generate.load_model`
334
335    *cutoff* is the polydispersity weight cutoff.
336    """
337    def __init__(self, data, model, cutoff=1e-5):
338        # type: (Data, KernelModel, float) -> None
339        self.model = model
340        self.cutoff = cutoff
341        # Note: _interpret_data defines the model attributes
342        self._interpret_data(data, model)
343
344    def __call__(self, **pars):
345        # type: (**float) -> np.ndarray
346        return self._calc_theory(pars, cutoff=self.cutoff)
347
348    def simulate_data(self, noise=None, **pars):
349        # type: (Optional[float], **float) -> None
350        """
351        Generate simulated data for the model.
352        """
353        Iq = self.__call__(**pars)
354        self._set_data(Iq, noise=noise)
355
356    def profile(self, **pars):
357        # type: (**float) -> None
358        """
359        Generate a plottable profile.
360        """
361        return call_profile(self.model.info, **pars)
362
363def main():
364    # type: () -> None
365    """
366    Program to evaluate a particular model at a set of q values.
367    """
368    import sys
369    from .data import empty_data1D, empty_data2D
370    from .core import load_model_info, build_model
371
372    if len(sys.argv) < 3:
373        print("usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ...")
374        sys.exit(1)
375    model_name = sys.argv[1]
376    call = sys.argv[2].upper()
377    if call != "ER_VR":
378        try:
379            values = [float(v) for v in call.split(',')]
380        except Exception:
381            values = []
382        if len(values) == 1:
383            q, = values
384            data = empty_data1D([q])
385        elif len(values) == 2:
386            qx, qy = values
387            data = empty_data2D([qx], [qy])
388        else:
389            print("use q or qx,qy or ER or VR")
390            sys.exit(1)
391    else:
392        data = empty_data1D([0.001])  # Data not used in ER/VR
393
394    model_info = load_model_info(model_name)
395    model = build_model(model_info)
396    calculator = DirectModel(data, model)
397    pars = dict((k, (float(v) if not k.endswith("_pd_type") else v))
398                for pair in sys.argv[3:]
399                for k, v in [pair.split('=')])
400    if call == "ER_VR":
401        ER = call_ER(model_info, pars)
402        VR = call_VR(model_info, pars)
403        print(ER, VR)
404    else:
405        Iq = calculator(**pars)
406        print(Iq[0])
407
408if __name__ == "__main__":
409    main()
Note: See TracBrowser for help on using the repository browser.