source: sasmodels/sasmodels/direct_model.py @ 304c775

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 304c775 was 304c775, checked in by Paul Kienzle <pkienzle@…>, 14 months ago

provide method for testing Fq results. Refs #1202.

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