[a84a0ca] | 1 | """ |
---|
| 2 | Wrap sasmodels for direct use by bumps. |
---|
| 3 | |
---|
| 4 | :class:`Model` is a wrapper for the sasmodels kernel which defines a |
---|
| 5 | bumps *Parameter* box for each kernel parameter. *Model* accepts keyword |
---|
| 6 | arguments to set the initial value for each parameter. |
---|
| 7 | |
---|
| 8 | :class:`Experiment` combines the *Model* function with a data file loaded by |
---|
| 9 | the sasview data loader. *Experiment* takes a *cutoff* parameter controlling |
---|
| 10 | how far the polydispersity integral extends. |
---|
| 11 | |
---|
| 12 | """ |
---|
[04dc697] | 13 | from __future__ import print_function |
---|
[a84a0ca] | 14 | |
---|
[40a87fa] | 15 | __all__ = ["Model", "Experiment"] |
---|
[a84a0ca] | 16 | |
---|
[7ae2b7f] | 17 | import numpy as np # type: ignore |
---|
[a84a0ca] | 18 | |
---|
| 19 | from .data import plot_theory |
---|
| 20 | from .direct_model import DataMixin |
---|
| 21 | |
---|
[04dc697] | 22 | try: |
---|
| 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] |
---|
| 28 | except ImportError: |
---|
| 29 | pass |
---|
| 30 | |
---|
| 31 | try: |
---|
| 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 |
---|
| 35 | except ImportError: |
---|
| 36 | pass |
---|
[a84a0ca] | 37 | |
---|
| 38 | |
---|
| 39 | def create_parameters(model_info, **kwargs): |
---|
[04dc697] | 40 | # type: (ModelInfo, **Union[float, str, Parameter]) -> Tuple[Dict[str, Parameter], Dict[str, str]] |
---|
[a84a0ca] | 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 |
---|
[04dc697] | 49 | value. The value can be a float, a bumps parameter, or in the case |
---|
| 50 | of the distribution type parameter, a string. |
---|
[a84a0ca] | 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 | """ |
---|
[04dc697] | 56 | pars = {} # type: Dict[str, Parameter] |
---|
| 57 | pd_types = {} # type: Dict[str, str] |
---|
[6d6508e] | 58 | for p in model_info.parameters.call_parameters: |
---|
[a84a0ca] | 59 | value = kwargs.pop(p.name, p.default) |
---|
| 60 | pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) |
---|
[d19962c] | 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) |
---|
[21b116f] | 70 | name = p.name + '_pd_type' |
---|
[04dc697] | 71 | pd_types[name] = str(kwargs.pop(name, 'gaussian')) |
---|
[a84a0ca] | 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 | |
---|
| 79 | class 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): |
---|
[04dc697] | 90 | # type: (KernelModel, **Dict[str, Union[float, Parameter]]) -> None |
---|
| 91 | self.sasmodel = model |
---|
[a84a0ca] | 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): |
---|
[04dc697] | 101 | # type: () -> Dict[str, Parameter] |
---|
[a84a0ca] | 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): |
---|
[04dc697] | 109 | # type: () -> Dict[str, Union[Parameter, str]] |
---|
[a84a0ca] | 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 | |
---|
| 118 | class 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 | """ |
---|
[04dc697] | 134 | _cache = None # type: Dict[str, np.ndarray] |
---|
[a84a0ca] | 135 | def __init__(self, data, model, cutoff=1e-5): |
---|
[04dc697] | 136 | # type: (Data, Model, float) -> None |
---|
[a84a0ca] | 137 | # remember inputs so we can inspect from outside |
---|
| 138 | self.model = model |
---|
| 139 | self.cutoff = cutoff |
---|
[04dc697] | 140 | self._interpret_data(data, model.sasmodel) |
---|
| 141 | self._cache = {} |
---|
[a84a0ca] | 142 | |
---|
| 143 | def update(self): |
---|
[04dc697] | 144 | # type: () -> None |
---|
[a84a0ca] | 145 | """ |
---|
| 146 | Call when model parameters have changed and theory needs to be |
---|
| 147 | recalculated. |
---|
| 148 | """ |
---|
[04dc697] | 149 | self._cache.clear() |
---|
[a84a0ca] | 150 | |
---|
| 151 | def numpoints(self): |
---|
[04dc697] | 152 | # type: () -> float |
---|
[a84a0ca] | 153 | """ |
---|
| 154 | Return the number of data points |
---|
| 155 | """ |
---|
| 156 | return len(self.Iq) |
---|
| 157 | |
---|
| 158 | def parameters(self): |
---|
[04dc697] | 159 | # type: () -> Dict[str, Parameter] |
---|
[a84a0ca] | 160 | """ |
---|
| 161 | Return a dictionary of parameters |
---|
| 162 | """ |
---|
| 163 | return self.model.parameters() |
---|
| 164 | |
---|
| 165 | def theory(self): |
---|
[04dc697] | 166 | # type: () -> np.ndarray |
---|
[a84a0ca] | 167 | """ |
---|
| 168 | Return the theory corresponding to the model parameters. |
---|
| 169 | |
---|
| 170 | This method uses lazy evaluation, and requires model.update() to be |
---|
| 171 | called when the parameters have changed. |
---|
| 172 | """ |
---|
| 173 | if 'theory' not in self._cache: |
---|
| 174 | pars = self.model.state() |
---|
| 175 | self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff) |
---|
| 176 | return self._cache['theory'] |
---|
| 177 | |
---|
| 178 | def residuals(self): |
---|
[04dc697] | 179 | # type: () -> np.ndarray |
---|
[a84a0ca] | 180 | """ |
---|
| 181 | Return theory minus data normalized by uncertainty. |
---|
| 182 | """ |
---|
| 183 | #if np.any(self.err ==0): print("zeros in err") |
---|
| 184 | return (self.theory() - self.Iq) / self.dIq |
---|
| 185 | |
---|
| 186 | def nllf(self): |
---|
[04dc697] | 187 | # type: () -> float |
---|
[a84a0ca] | 188 | """ |
---|
| 189 | Return the negative log likelihood of seeing data given the model |
---|
| 190 | parameters, up to a normalizing constant which depends on the data |
---|
| 191 | uncertainty. |
---|
| 192 | """ |
---|
| 193 | delta = self.residuals() |
---|
| 194 | #if np.any(np.isnan(R)): print("NaN in residuals") |
---|
[04dc697] | 195 | return 0.5 * np.sum(delta**2) |
---|
[a84a0ca] | 196 | |
---|
| 197 | #def __call__(self): |
---|
| 198 | # return 2 * self.nllf() / self.dof |
---|
| 199 | |
---|
| 200 | def plot(self, view='log'): |
---|
[04dc697] | 201 | # type: (str) -> None |
---|
[a84a0ca] | 202 | """ |
---|
| 203 | Plot the data and residuals. |
---|
| 204 | """ |
---|
| 205 | data, theory, resid = self._data, self.theory(), self.residuals() |
---|
[40a87fa] | 206 | plot_theory(data, theory, resid, view, Iq_calc=self.Iq_calc) |
---|
[a84a0ca] | 207 | |
---|
| 208 | def simulate_data(self, noise=None): |
---|
[04dc697] | 209 | # type: (float) -> None |
---|
[a84a0ca] | 210 | """ |
---|
| 211 | Generate simulated data. |
---|
| 212 | """ |
---|
| 213 | Iq = self.theory() |
---|
| 214 | self._set_data(Iq, noise) |
---|
| 215 | |
---|
| 216 | def save(self, basename): |
---|
[04dc697] | 217 | # type: (str) -> None |
---|
[a84a0ca] | 218 | """ |
---|
| 219 | Save the model parameters and data into a file. |
---|
| 220 | |
---|
| 221 | Not Implemented. |
---|
| 222 | """ |
---|
[9217ef8] | 223 | if self.data_type == "sesans": |
---|
| 224 | np.savetxt(basename+".dat", np.array([self._data.x, self.theory()]).T) |
---|
[a84a0ca] | 225 | |
---|
| 226 | def __getstate__(self): |
---|
[04dc697] | 227 | # type: () -> Dict[str, Any] |
---|
[a84a0ca] | 228 | # Can't pickle gpu functions, so instead make them lazy |
---|
| 229 | state = self.__dict__.copy() |
---|
| 230 | state['_kernel'] = None |
---|
| 231 | return state |
---|
| 232 | |
---|
| 233 | def __setstate__(self, state): |
---|
[04dc697] | 234 | # type: (Dict[str, Any]) -> None |
---|
[a84a0ca] | 235 | # pylint: disable=attribute-defined-outside-init |
---|
| 236 | self.__dict__ = state |
---|
[04dc697] | 237 | |
---|