[ff7119b] | 1 | """ |
---|
[aa4946b] | 2 | Wrap sasmodels for direct use by bumps. |
---|
[346bc88] | 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 | |
---|
[37a7252] | 8 | :class:`Experiment` combines the *Model* function with a data file loaded by |
---|
| 9 | the sasview data loader. *Experiment* takes a *cutoff* parameter controlling |
---|
[346bc88] | 10 | how far the polydispersity integral extends. |
---|
| 11 | |
---|
[ff7119b] | 12 | """ |
---|
[aa4946b] | 13 | |
---|
[346bc88] | 14 | import warnings |
---|
[14de349] | 15 | |
---|
[aa4946b] | 16 | import numpy as np |
---|
| 17 | |
---|
[7cf2cfd] | 18 | from .data import plot_theory |
---|
| 19 | from .direct_model import DataMixin |
---|
[14de349] | 20 | |
---|
[190fc2b] | 21 | __all__ = [ |
---|
| 22 | "Model", "Experiment", |
---|
| 23 | ] |
---|
| 24 | |
---|
[346bc88] | 25 | # CRUFT: old style bumps wrapper which doesn't separate data and model |
---|
[37a7252] | 26 | # pylint: disable=invalid-name |
---|
[346bc88] | 27 | def BumpsModel(data, model, cutoff=1e-5, **kw): |
---|
[37a7252] | 28 | r""" |
---|
| 29 | Bind a model to data, along with a polydispersity cutoff. |
---|
| 30 | |
---|
| 31 | *data* is a :class:`data.Data1D`, :class:`data.Data2D` or |
---|
| 32 | :class:`data.Sesans` object. Use :func:`data.empty_data1D` or |
---|
| 33 | :func:`data.empty_data2D` to define $q, \Delta q$ calculation |
---|
| 34 | points for displaying the SANS curve when there is no measured data. |
---|
| 35 | |
---|
| 36 | *model* is a runnable module as returned from :func:`core.load_model`. |
---|
| 37 | |
---|
| 38 | *cutoff* is the polydispersity weight cutoff. |
---|
| 39 | |
---|
| 40 | Any additional *key=value* pairs are model dependent parameters. |
---|
| 41 | |
---|
| 42 | Returns an :class:`Experiment` object. |
---|
| 43 | |
---|
| 44 | Note that the usual Bumps semantics is not fully supported, since |
---|
| 45 | assigning *M.name = parameter* on the returned experiment object |
---|
| 46 | does not set that parameter in the model. Range setting will still |
---|
| 47 | work as expected though. |
---|
| 48 | |
---|
| 49 | .. deprecated:: 0.1 |
---|
| 50 | Use :class:`Experiment` instead. |
---|
| 51 | """ |
---|
[346bc88] | 52 | warnings.warn("Use of BumpsModel is deprecated. Use bumps_model.Experiment instead.") |
---|
[37a7252] | 53 | |
---|
| 54 | # Create the model and experiment |
---|
[346bc88] | 55 | model = Model(model, **kw) |
---|
| 56 | experiment = Experiment(data=data, model=model, cutoff=cutoff) |
---|
[37a7252] | 57 | |
---|
| 58 | # Copy the model parameters up to the experiment object. |
---|
| 59 | for k, v in model.parameters().items(): |
---|
| 60 | setattr(experiment, k, v) |
---|
[346bc88] | 61 | return experiment |
---|
| 62 | |
---|
[37a7252] | 63 | |
---|
[ec7e360] | 64 | def create_parameters(model_info, **kwargs): |
---|
[37a7252] | 65 | """ |
---|
| 66 | Generate Bumps parameters from the model info. |
---|
| 67 | |
---|
| 68 | *model_info* is returned from :func:`generate.model_info` on the |
---|
| 69 | model definition module. |
---|
| 70 | |
---|
| 71 | Any additional *key=value* pairs are initial values for the parameters |
---|
| 72 | to the models. Uninitialized parameters will use the model default |
---|
| 73 | value. |
---|
| 74 | |
---|
| 75 | Returns a dictionary of *{name: Parameter}* containing the bumps |
---|
| 76 | parameters for each model parameter, and a dictionary of |
---|
| 77 | *{name: str}* containing the polydispersity distribution types. |
---|
| 78 | """ |
---|
[ec7e360] | 79 | # lazy import; this allows the doc builder and nosetests to run even |
---|
| 80 | # when bumps is not on the path. |
---|
| 81 | from bumps.names import Parameter |
---|
| 82 | |
---|
| 83 | pars = {} |
---|
| 84 | for p in model_info['parameters']: |
---|
[fcd7bbd] | 85 | value = kwargs.pop([p.name, p.default]) |
---|
| 86 | pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) |
---|
[37a7252] | 87 | for name in model_info['partype']['pd-2d']: |
---|
[ec7e360] | 88 | for xpart, xdefault, xlimits in [ |
---|
[fcd7bbd] | 89 | ('_pd', 0., pars[name].limits), |
---|
[37a7252] | 90 | ('_pd_n', 35., (0, 1000)), |
---|
| 91 | ('_pd_nsigma', 3., (0, 10)), |
---|
| 92 | ]: |
---|
[ec7e360] | 93 | xname = name + xpart |
---|
| 94 | xvalue = kwargs.pop(xname, xdefault) |
---|
| 95 | pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) |
---|
| 96 | |
---|
| 97 | pd_types = {} |
---|
[37a7252] | 98 | for name in model_info['partype']['pd-2d']: |
---|
[ec7e360] | 99 | xname = name + '_pd_type' |
---|
| 100 | xvalue = kwargs.pop(xname, 'gaussian') |
---|
| 101 | pd_types[xname] = xvalue |
---|
| 102 | |
---|
| 103 | if kwargs: # args not corresponding to parameters |
---|
| 104 | raise TypeError("unexpected parameters: %s" |
---|
| 105 | % (", ".join(sorted(kwargs.keys())))) |
---|
| 106 | |
---|
| 107 | return pars, pd_types |
---|
[346bc88] | 108 | |
---|
| 109 | class Model(object): |
---|
[37a7252] | 110 | """ |
---|
| 111 | Bumps wrapper for a SAS model. |
---|
| 112 | |
---|
| 113 | *model* is a runnable module as returned from :func:`core.load_model`. |
---|
| 114 | |
---|
| 115 | *cutoff* is the polydispersity weight cutoff. |
---|
| 116 | |
---|
| 117 | Any additional *key=value* pairs are model dependent parameters. |
---|
| 118 | """ |
---|
[ec7e360] | 119 | def __init__(self, model, **kwargs): |
---|
[7cf2cfd] | 120 | self._sasmodel = model |
---|
[ec7e360] | 121 | pars, pd_types = create_parameters(model.info, **kwargs) |
---|
[37a7252] | 122 | for k, v in pars.items(): |
---|
[ec7e360] | 123 | setattr(self, k, v) |
---|
[37a7252] | 124 | for k, v in pd_types.items(): |
---|
[ec7e360] | 125 | setattr(self, k, v) |
---|
| 126 | self._parameter_names = list(pars.keys()) |
---|
| 127 | self._pd_type_names = list(pd_types.keys()) |
---|
[346bc88] | 128 | |
---|
| 129 | def parameters(self): |
---|
| 130 | """ |
---|
[37a7252] | 131 | Return a dictionary of parameters objects for the parameters, |
---|
| 132 | excluding polydispersity distribution type. |
---|
[346bc88] | 133 | """ |
---|
| 134 | return dict((k, getattr(self, k)) for k in self._parameter_names) |
---|
| 135 | |
---|
[ec7e360] | 136 | def state(self): |
---|
[37a7252] | 137 | """ |
---|
| 138 | Return a dictionary of current values for all the parameters, |
---|
| 139 | including polydispersity distribution type. |
---|
| 140 | """ |
---|
[ec7e360] | 141 | pars = dict((k, getattr(self, k).value) for k in self._parameter_names) |
---|
| 142 | pars.update((k, getattr(self, k)) for k in self._pd_type_names) |
---|
| 143 | return pars |
---|
[7cf2cfd] | 144 | |
---|
| 145 | class Experiment(DataMixin): |
---|
[37a7252] | 146 | r""" |
---|
| 147 | Bumps wrapper for a SAS experiment. |
---|
[ff7119b] | 148 | |
---|
[37a7252] | 149 | *data* is a :class:`data.Data1D`, :class:`data.Data2D` or |
---|
| 150 | :class:`data.Sesans` object. Use :func:`data.empty_data1D` or |
---|
| 151 | :func:`data.empty_data2D` to define $q, \Delta q$ calculation |
---|
| 152 | points for displaying the SANS curve when there is no measured data. |
---|
[ff7119b] | 153 | |
---|
[37a7252] | 154 | *model* is a :class:`Model` object. |
---|
[ff7119b] | 155 | |
---|
| 156 | *cutoff* is the integration cutoff, which avoids computing the |
---|
| 157 | the SAS model where the polydispersity weight is low. |
---|
| 158 | |
---|
[37a7252] | 159 | The resulting model can be used directly in a Bumps FitProblem call. |
---|
[ff7119b] | 160 | """ |
---|
[346bc88] | 161 | def __init__(self, data, model, cutoff=1e-5): |
---|
[14de349] | 162 | |
---|
[87985ca] | 163 | # remember inputs so we can inspect from outside |
---|
| 164 | self.model = model |
---|
[abb22f4] | 165 | self.cutoff = cutoff |
---|
[7cf2cfd] | 166 | self._interpret_data(data, model._sasmodel) |
---|
[14de349] | 167 | self.update() |
---|
| 168 | |
---|
| 169 | def update(self): |
---|
[37a7252] | 170 | """ |
---|
| 171 | Call when model parameters have changed and theory needs to be |
---|
| 172 | recalculated. |
---|
| 173 | """ |
---|
[14de349] | 174 | self._cache = {} |
---|
| 175 | |
---|
| 176 | def numpoints(self): |
---|
[7e224c2] | 177 | """ |
---|
[37a7252] | 178 | Return the number of data points |
---|
[7e224c2] | 179 | """ |
---|
[3c56da87] | 180 | return len(self.Iq) |
---|
[14de349] | 181 | |
---|
| 182 | def parameters(self): |
---|
[7e224c2] | 183 | """ |
---|
[346bc88] | 184 | Return a dictionary of parameters |
---|
[7e224c2] | 185 | """ |
---|
[346bc88] | 186 | return self.model.parameters() |
---|
[14de349] | 187 | |
---|
| 188 | def theory(self): |
---|
[37a7252] | 189 | """ |
---|
| 190 | Return the theory corresponding to the model parameters. |
---|
| 191 | |
---|
| 192 | This method uses lazy evaluation, and requires model.update() to be |
---|
| 193 | called when the parameters have changed. |
---|
| 194 | """ |
---|
[14de349] | 195 | if 'theory' not in self._cache: |
---|
[ec7e360] | 196 | pars = self.model.state() |
---|
[7cf2cfd] | 197 | self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff) |
---|
[14de349] | 198 | return self._cache['theory'] |
---|
| 199 | |
---|
| 200 | def residuals(self): |
---|
[37a7252] | 201 | """ |
---|
| 202 | Return theory minus data normalized by uncertainty. |
---|
| 203 | """ |
---|
[9404dd3] | 204 | #if np.any(self.err ==0): print("zeros in err") |
---|
[346bc88] | 205 | return (self.theory() - self.Iq) / self.dIq |
---|
[14de349] | 206 | |
---|
| 207 | def nllf(self): |
---|
[37a7252] | 208 | """ |
---|
| 209 | Return the negative log likelihood of seeing data given the model |
---|
| 210 | parameters, up to a normalizing constant which depends on the data |
---|
| 211 | uncertainty. |
---|
| 212 | """ |
---|
[3c56da87] | 213 | delta = self.residuals() |
---|
[9404dd3] | 214 | #if np.any(np.isnan(R)): print("NaN in residuals") |
---|
[3c56da87] | 215 | return 0.5 * np.sum(delta ** 2) |
---|
[14de349] | 216 | |
---|
[3c56da87] | 217 | #def __call__(self): |
---|
| 218 | # return 2 * self.nllf() / self.dof |
---|
[14de349] | 219 | |
---|
| 220 | def plot(self, view='log'): |
---|
[c97724e] | 221 | """ |
---|
| 222 | Plot the data and residuals. |
---|
| 223 | """ |
---|
[7cf2cfd] | 224 | data, theory, resid = self._data, self.theory(), self.residuals() |
---|
| 225 | plot_theory(data, theory, resid, view) |
---|
[c97724e] | 226 | |
---|
| 227 | def simulate_data(self, noise=None): |
---|
[37a7252] | 228 | """ |
---|
| 229 | Generate simulated data. |
---|
| 230 | """ |
---|
[7cf2cfd] | 231 | Iq = self.theory() |
---|
| 232 | self._set_data(Iq, noise) |
---|
[14de349] | 233 | |
---|
| 234 | def save(self, basename): |
---|
[37a7252] | 235 | """ |
---|
| 236 | Save the model parameters and data into a file. |
---|
| 237 | |
---|
| 238 | Not Implemented. |
---|
| 239 | """ |
---|
[14de349] | 240 | pass |
---|
| 241 | |
---|
[abb22f4] | 242 | def __getstate__(self): |
---|
| 243 | # Can't pickle gpu functions, so instead make them lazy |
---|
| 244 | state = self.__dict__.copy() |
---|
[7cf2cfd] | 245 | state['_kernel'] = None |
---|
[abb22f4] | 246 | return state |
---|
| 247 | |
---|
| 248 | def __setstate__(self, state): |
---|
[3c56da87] | 249 | # pylint: disable=attribute-defined-outside-init |
---|
[abb22f4] | 250 | self.__dict__ = state |
---|