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 | """ |
---|
13 | |
---|
14 | __all__ = [ |
---|
15 | "Model", "Experiment", |
---|
16 | ] |
---|
17 | |
---|
18 | import warnings |
---|
19 | |
---|
20 | import numpy as np |
---|
21 | |
---|
22 | from .data import plot_theory |
---|
23 | from .direct_model import DataMixin |
---|
24 | |
---|
25 | # CRUFT: old style bumps wrapper which doesn't separate data and model |
---|
26 | # pylint: disable=invalid-name |
---|
27 | def BumpsModel(data, model, cutoff=1e-5, **kw): |
---|
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 | """ |
---|
52 | warnings.warn("Use of BumpsModel is deprecated. Use bumps_model.Experiment instead.") |
---|
53 | |
---|
54 | # Create the model and experiment |
---|
55 | model = Model(model, **kw) |
---|
56 | experiment = Experiment(data=data, model=model, cutoff=cutoff) |
---|
57 | |
---|
58 | # Copy the model parameters up to the experiment object. |
---|
59 | for k, v in model.parameters().items(): |
---|
60 | setattr(experiment, k, v) |
---|
61 | return experiment |
---|
62 | |
---|
63 | |
---|
64 | def create_parameters(model_info, **kwargs): |
---|
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 | """ |
---|
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']: |
---|
85 | name, default, limits = p[0], p[2], p[3] |
---|
86 | value = kwargs.pop(name, default) |
---|
87 | pars[name] = Parameter.default(value, name=name, limits=limits) |
---|
88 | for name in model_info['partype']['pd-2d']: |
---|
89 | for xpart, xdefault, xlimits in [ |
---|
90 | ('_pd', 0., limits), |
---|
91 | ('_pd_n', 35., (0, 1000)), |
---|
92 | ('_pd_nsigma', 3., (0, 10)), |
---|
93 | ]: |
---|
94 | xname = name + xpart |
---|
95 | xvalue = kwargs.pop(xname, xdefault) |
---|
96 | pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) |
---|
97 | |
---|
98 | pd_types = {} |
---|
99 | for name in model_info['partype']['pd-2d']: |
---|
100 | xname = name + '_pd_type' |
---|
101 | xvalue = kwargs.pop(xname, 'gaussian') |
---|
102 | pd_types[xname] = xvalue |
---|
103 | |
---|
104 | if kwargs: # args not corresponding to parameters |
---|
105 | raise TypeError("unexpected parameters: %s" |
---|
106 | % (", ".join(sorted(kwargs.keys())))) |
---|
107 | |
---|
108 | return pars, pd_types |
---|
109 | |
---|
110 | class Model(object): |
---|
111 | """ |
---|
112 | Bumps wrapper for a SAS model. |
---|
113 | |
---|
114 | *model* is a runnable module as returned from :func:`core.load_model`. |
---|
115 | |
---|
116 | *cutoff* is the polydispersity weight cutoff. |
---|
117 | |
---|
118 | Any additional *key=value* pairs are model dependent parameters. |
---|
119 | """ |
---|
120 | def __init__(self, model, **kwargs): |
---|
121 | self._sasmodel = model |
---|
122 | pars, pd_types = create_parameters(model.info, **kwargs) |
---|
123 | for k, v in pars.items(): |
---|
124 | setattr(self, k, v) |
---|
125 | for k, v in pd_types.items(): |
---|
126 | setattr(self, k, v) |
---|
127 | self._parameter_names = list(pars.keys()) |
---|
128 | self._pd_type_names = list(pd_types.keys()) |
---|
129 | |
---|
130 | def parameters(self): |
---|
131 | """ |
---|
132 | Return a dictionary of parameters objects for the parameters, |
---|
133 | excluding polydispersity distribution type. |
---|
134 | """ |
---|
135 | return dict((k, getattr(self, k)) for k in self._parameter_names) |
---|
136 | |
---|
137 | def state(self): |
---|
138 | """ |
---|
139 | Return a dictionary of current values for all the parameters, |
---|
140 | including polydispersity distribution type. |
---|
141 | """ |
---|
142 | pars = dict((k, getattr(self, k).value) for k in self._parameter_names) |
---|
143 | pars.update((k, getattr(self, k)) for k in self._pd_type_names) |
---|
144 | return pars |
---|
145 | |
---|
146 | class Experiment(DataMixin): |
---|
147 | r""" |
---|
148 | Bumps wrapper for a SAS experiment. |
---|
149 | |
---|
150 | *data* is a :class:`data.Data1D`, :class:`data.Data2D` or |
---|
151 | :class:`data.Sesans` object. Use :func:`data.empty_data1D` or |
---|
152 | :func:`data.empty_data2D` to define $q, \Delta q$ calculation |
---|
153 | points for displaying the SANS curve when there is no measured data. |
---|
154 | |
---|
155 | *model* is a :class:`Model` object. |
---|
156 | |
---|
157 | *cutoff* is the integration cutoff, which avoids computing the |
---|
158 | the SAS model where the polydispersity weight is low. |
---|
159 | |
---|
160 | The resulting model can be used directly in a Bumps FitProblem call. |
---|
161 | """ |
---|
162 | def __init__(self, data, model, cutoff=1e-5): |
---|
163 | |
---|
164 | # remember inputs so we can inspect from outside |
---|
165 | self.model = model |
---|
166 | self.cutoff = cutoff |
---|
167 | self._interpret_data(data, model._sasmodel) |
---|
168 | self.update() |
---|
169 | |
---|
170 | def update(self): |
---|
171 | """ |
---|
172 | Call when model parameters have changed and theory needs to be |
---|
173 | recalculated. |
---|
174 | """ |
---|
175 | self._cache = {} |
---|
176 | |
---|
177 | def numpoints(self): |
---|
178 | """ |
---|
179 | Return the number of data points |
---|
180 | """ |
---|
181 | return len(self.Iq) |
---|
182 | |
---|
183 | def parameters(self): |
---|
184 | """ |
---|
185 | Return a dictionary of parameters |
---|
186 | """ |
---|
187 | return self.model.parameters() |
---|
188 | |
---|
189 | def theory(self): |
---|
190 | """ |
---|
191 | Return the theory corresponding to the model parameters. |
---|
192 | |
---|
193 | This method uses lazy evaluation, and requires model.update() to be |
---|
194 | called when the parameters have changed. |
---|
195 | """ |
---|
196 | if 'theory' not in self._cache: |
---|
197 | pars = self.model.state() |
---|
198 | self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff) |
---|
199 | return self._cache['theory'] |
---|
200 | |
---|
201 | def residuals(self): |
---|
202 | """ |
---|
203 | Return theory minus data normalized by uncertainty. |
---|
204 | """ |
---|
205 | #if np.any(self.err ==0): print("zeros in err") |
---|
206 | return (self.theory() - self.Iq) / self.dIq |
---|
207 | |
---|
208 | def nllf(self): |
---|
209 | """ |
---|
210 | Return the negative log likelihood of seeing data given the model |
---|
211 | parameters, up to a normalizing constant which depends on the data |
---|
212 | uncertainty. |
---|
213 | """ |
---|
214 | delta = self.residuals() |
---|
215 | #if np.any(np.isnan(R)): print("NaN in residuals") |
---|
216 | return 0.5 * np.sum(delta ** 2) |
---|
217 | |
---|
218 | #def __call__(self): |
---|
219 | # return 2 * self.nllf() / self.dof |
---|
220 | |
---|
221 | def plot(self, view='log'): |
---|
222 | """ |
---|
223 | Plot the data and residuals. |
---|
224 | """ |
---|
225 | data, theory, resid = self._data, self.theory(), self.residuals() |
---|
226 | plot_theory(data, theory, resid, view) |
---|
227 | |
---|
228 | def simulate_data(self, noise=None): |
---|
229 | """ |
---|
230 | Generate simulated data. |
---|
231 | """ |
---|
232 | Iq = self.theory() |
---|
233 | self._set_data(Iq, noise) |
---|
234 | |
---|
235 | def save(self, basename): |
---|
236 | """ |
---|
237 | Save the model parameters and data into a file. |
---|
238 | |
---|
239 | Not Implemented. |
---|
240 | """ |
---|
241 | pass |
---|
242 | |
---|
243 | def __getstate__(self): |
---|
244 | # Can't pickle gpu functions, so instead make them lazy |
---|
245 | state = self.__dict__.copy() |
---|
246 | state['_kernel'] = None |
---|
247 | return state |
---|
248 | |
---|
249 | def __setstate__(self, state): |
---|
250 | # pylint: disable=attribute-defined-outside-init |
---|
251 | self.__dict__ = state |
---|