source: sasmodels/sasmodels/bumps_model.py @ 74b0495

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

update oriented sesans and simultaneous fit examples

  • Property mode set to 100644
File size: 8.0 KB
Line 
1"""
2Wrap sasmodels for direct use by bumps.
3
4:class:`Model` is a wrapper for the sasmodels kernel which defines a
5bumps *Parameter* box for each kernel parameter.  *Model* accepts keyword
6arguments to set the initial value for each parameter.
7
8:class:`Experiment` combines the *Model* function with a data file loaded by
9the sasview data loader.  *Experiment* takes a *cutoff* parameter controlling
10how far the polydispersity integral extends.
11
12"""
13from __future__ import print_function
14
15__all__ = ["Model", "Experiment"]
16
17import numpy as np  # type: ignore
18
19from .data import plot_theory
20from .direct_model import DataMixin
21
22try:
23    from typing import Dict, Union, Tuple, Any
24    from .data import Data1D, Data2D
25    from .kernel import KernelModel
26    from .modelinfo import ModelInfo
27    Data = Union[Data1D, Data2D]
28except ImportError:
29    pass
30
31try:
32    # Optional import. This allows the doc builder and nosetests to run even
33    # when bumps is not on the path.
34    from bumps.names import Parameter # type: ignore
35except ImportError:
36    pass
37
38
39def create_parameters(model_info, **kwargs):
40    # type: (ModelInfo, **Union[float, str, Parameter]) -> Tuple[Dict[str, Parameter], Dict[str, str]]
41    """
42    Generate Bumps parameters from the model info.
43
44    *model_info* is returned from :func:`generate.model_info` on the
45    model definition module.
46
47    Any additional *key=value* pairs are initial values for the parameters
48    to the models.  Uninitialized parameters will use the model default
49    value.  The value can be a float, a bumps parameter, or in the case
50    of the distribution type parameter, a string.
51
52    Returns a dictionary of *{name: Parameter}* containing the bumps
53    parameters for each model parameter, and a dictionary of
54    *{name: str}* containing the polydispersity distribution types.
55    """
56    pars = {}     # type: Dict[str, Parameter]
57    pd_types = {} # type: Dict[str, str]
58    for p in model_info.parameters.call_parameters:
59        value = kwargs.pop(p.name, p.default)
60        pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits)
61        if p.polydisperse:
62            for part, default, limits in [
63                    ('_pd', 0., pars[p.name].limits),
64                    ('_pd_n', 35., (0, 1000)),
65                    ('_pd_nsigma', 3., (0, 10)),
66                ]:
67                name = p.name + part
68                value = kwargs.pop(name, default)
69                pars[name] = Parameter.default(value, name=name, limits=limits)
70            name = p.name + '_pd_type'
71            pd_types[name] = str(kwargs.pop(name, 'gaussian'))
72
73    if kwargs:  # args not corresponding to parameters
74        raise TypeError("unexpected parameters: %s"
75                        % (", ".join(sorted(kwargs.keys()))))
76
77    return pars, pd_types
78
79class Model(object):
80    """
81    Bumps wrapper for a SAS model.
82
83    *model* is a runnable module as returned from :func:`core.load_model`.
84
85    *cutoff* is the polydispersity weight cutoff.
86
87    Any additional *key=value* pairs are model dependent parameters.
88    """
89    def __init__(self, model, **kwargs):
90        # type: (KernelModel, **Dict[str, Union[float, Parameter]]) -> None
91        self.sasmodel = model
92        pars, pd_types = create_parameters(model.info, **kwargs)
93        for k, v in pars.items():
94            setattr(self, k, v)
95        for k, v in pd_types.items():
96            setattr(self, k, v)
97        self._parameter_names = list(pars.keys())
98        self._pd_type_names = list(pd_types.keys())
99
100    def parameters(self):
101        # type: () -> Dict[str, Parameter]
102        """
103        Return a dictionary of parameters objects for the parameters,
104        excluding polydispersity distribution type.
105        """
106        return dict((k, getattr(self, k)) for k in self._parameter_names)
107
108    def state(self):
109        # type: () -> Dict[str, Union[Parameter, str]]
110        """
111        Return a dictionary of current values for all the parameters,
112        including polydispersity distribution type.
113        """
114        pars = dict((k, getattr(self, k).value) for k in self._parameter_names)
115        pars.update((k, getattr(self, k)) for k in self._pd_type_names)
116        return pars
117
118class Experiment(DataMixin):
119    r"""
120    Bumps wrapper for a SAS experiment.
121
122    *data* is a :class:`data.Data1D`, :class:`data.Data2D` or
123    :class:`data.Sesans` object.  Use :func:`data.empty_data1D` or
124    :func:`data.empty_data2D` to define $q, \Delta q$ calculation
125    points for displaying the SANS curve when there is no measured data.
126
127    *model* is a :class:`Model` object.
128
129    *cutoff* is the integration cutoff, which avoids computing the
130    the SAS model where the polydispersity weight is low.
131
132    The resulting model can be used directly in a Bumps FitProblem call.
133    """
134    _cache = None # type: Dict[str, np.ndarray]
135    def __init__(self, data, model, cutoff=1e-5, name=None):
136        # type: (Data, Model, float) -> None
137        # remember inputs so we can inspect from outside
138        self.name = data.filename if name is None else name
139        self.model = model
140        self.cutoff = cutoff
141        self._interpret_data(data, model.sasmodel)
142        self._cache = {}
143
144    def update(self):
145        # type: () -> None
146        """
147        Call when model parameters have changed and theory needs to be
148        recalculated.
149        """
150        self._cache.clear()
151
152    def numpoints(self):
153        # type: () -> float
154        """
155        Return the number of data points
156        """
157        return len(self.Iq)
158
159    def parameters(self):
160        # type: () -> Dict[str, Parameter]
161        """
162        Return a dictionary of parameters
163        """
164        return self.model.parameters()
165
166    def theory(self):
167        # type: () -> np.ndarray
168        """
169        Return the theory corresponding to the model parameters.
170
171        This method uses lazy evaluation, and requires model.update() to be
172        called when the parameters have changed.
173        """
174        if 'theory' not in self._cache:
175            pars = self.model.state()
176            self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff)
177        return self._cache['theory']
178
179    def residuals(self):
180        # type: () -> np.ndarray
181        """
182        Return theory minus data normalized by uncertainty.
183        """
184        #if np.any(self.err ==0): print("zeros in err")
185        return (self.theory() - self.Iq) / self.dIq
186
187    def nllf(self):
188        # type: () -> float
189        """
190        Return the negative log likelihood of seeing data given the model
191        parameters, up to a normalizing constant which depends on the data
192        uncertainty.
193        """
194        delta = self.residuals()
195        #if np.any(np.isnan(R)): print("NaN in residuals")
196        return 0.5 * np.sum(delta**2)
197
198    #def __call__(self):
199    #    return 2 * self.nllf() / self.dof
200
201    def plot(self, view='log'):
202        # type: (str) -> None
203        """
204        Plot the data and residuals.
205        """
206        data, theory, resid = self._data, self.theory(), self.residuals()
207        # TODO: hack to display oriented usans 2-D pattern
208        Iq_calc = self.Iq_calc if isinstance(self.Iq_calc, tuple) else None
209        plot_theory(data, theory, resid, view, Iq_calc=Iq_calc)
210
211    def simulate_data(self, noise=None):
212        # type: (float) -> None
213        """
214        Generate simulated data.
215        """
216        Iq = self.theory()
217        self._set_data(Iq, noise)
218
219    def save(self, basename):
220        # type: (str) -> None
221        """
222        Save the model parameters and data into a file.
223
224        Not Implemented.
225        """
226        if self.data_type == "sesans":
227            np.savetxt(basename+".dat", np.array([self._data.x, self.theory()]).T)
228
229    def __getstate__(self):
230        # type: () -> Dict[str, Any]
231        # Can't pickle gpu functions, so instead make them lazy
232        state = self.__dict__.copy()
233        state['_kernel'] = None
234        return state
235
236    def __setstate__(self, state):
237        # type: (Dict[str, Any]) -> None
238        # pylint: disable=attribute-defined-outside-init
239        self.__dict__ = state
240
Note: See TracBrowser for help on using the repository browser.