[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 | |
---|
| 8 | :class:`Experiment` combines the *Model* function with a data file loaded by the |
---|
| 9 | sasview data loader. *Experiment* takes a *cutoff* parameter controlling |
---|
| 10 | how far the polydispersity integral extends. |
---|
| 11 | |
---|
[ff7119b] | 12 | """ |
---|
[aa4946b] | 13 | |
---|
[14de349] | 14 | import datetime |
---|
[346bc88] | 15 | import warnings |
---|
[14de349] | 16 | |
---|
[aa4946b] | 17 | import numpy as np |
---|
| 18 | |
---|
[7cf2cfd] | 19 | from . import sesans |
---|
| 20 | from . import weights |
---|
| 21 | from .data import plot_theory |
---|
| 22 | from .direct_model import DataMixin |
---|
[14de349] | 23 | |
---|
[346bc88] | 24 | # CRUFT: old style bumps wrapper which doesn't separate data and model |
---|
| 25 | def BumpsModel(data, model, cutoff=1e-5, **kw): |
---|
| 26 | warnings.warn("Use of BumpsModel is deprecated. Use bumps_model.Experiment instead.") |
---|
| 27 | model = Model(model, **kw) |
---|
| 28 | experiment = Experiment(data=data, model=model, cutoff=cutoff) |
---|
| 29 | for k in model._parameter_names: |
---|
| 30 | setattr(experiment, k, getattr(model, k)) |
---|
| 31 | return experiment |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | class Model(object): |
---|
[7cf2cfd] | 35 | def __init__(self, model, **kw): |
---|
[0e9048f] | 36 | # lazy import; this allows the doc builder and nosetests to run even |
---|
| 37 | # when bumps is not on the path. |
---|
| 38 | from bumps.names import Parameter |
---|
| 39 | |
---|
[7cf2cfd] | 40 | self._sasmodel = model |
---|
| 41 | partype = model.info['partype'] |
---|
[346bc88] | 42 | |
---|
| 43 | pars = [] |
---|
[7cf2cfd] | 44 | for p in model.info['parameters']: |
---|
[346bc88] | 45 | name, default, limits = p[0], p[2], p[3] |
---|
| 46 | value = kw.pop(name, default) |
---|
| 47 | setattr(self, name, Parameter.default(value, name=name, limits=limits)) |
---|
| 48 | pars.append(name) |
---|
| 49 | for name in partype['pd-2d']: |
---|
| 50 | for xpart, xdefault, xlimits in [ |
---|
| 51 | ('_pd', 0, limits), |
---|
| 52 | ('_pd_n', 35, (0, 1000)), |
---|
| 53 | ('_pd_nsigma', 3, (0, 10)), |
---|
| 54 | ('_pd_type', 'gaussian', None), |
---|
| 55 | ]: |
---|
| 56 | xname = name + xpart |
---|
| 57 | xvalue = kw.pop(xname, xdefault) |
---|
| 58 | if xlimits is not None: |
---|
| 59 | xvalue = Parameter.default(xvalue, name=xname, limits=xlimits) |
---|
| 60 | pars.append(xname) |
---|
| 61 | setattr(self, xname, xvalue) |
---|
| 62 | self._parameter_names = pars |
---|
| 63 | if kw: |
---|
| 64 | raise TypeError("unexpected parameters: %s" |
---|
| 65 | % (", ".join(sorted(kw.keys())))) |
---|
| 66 | |
---|
| 67 | def parameters(self): |
---|
| 68 | """ |
---|
| 69 | Return a dictionary of parameters |
---|
| 70 | """ |
---|
| 71 | return dict((k, getattr(self, k)) for k in self._parameter_names) |
---|
| 72 | |
---|
[7cf2cfd] | 73 | |
---|
| 74 | class Experiment(DataMixin): |
---|
[ff7119b] | 75 | """ |
---|
| 76 | Return a bumps wrapper for a SAS model. |
---|
| 77 | |
---|
| 78 | *data* is the data to be fitted. |
---|
| 79 | |
---|
[aa4946b] | 80 | *model* is the SAS model from :func:`core.load_model`. |
---|
[ff7119b] | 81 | |
---|
| 82 | *cutoff* is the integration cutoff, which avoids computing the |
---|
| 83 | the SAS model where the polydispersity weight is low. |
---|
| 84 | |
---|
| 85 | Model parameters can be initialized with additional keyword |
---|
| 86 | arguments, or by assigning to model.parameter_name.value. |
---|
| 87 | |
---|
| 88 | The resulting bumps model can be used directly in a FitProblem call. |
---|
| 89 | """ |
---|
[346bc88] | 90 | def __init__(self, data, model, cutoff=1e-5): |
---|
[14de349] | 91 | |
---|
[87985ca] | 92 | # remember inputs so we can inspect from outside |
---|
| 93 | self.model = model |
---|
[abb22f4] | 94 | self.cutoff = cutoff |
---|
[7cf2cfd] | 95 | self._interpret_data(data, model._sasmodel) |
---|
[14de349] | 96 | self.update() |
---|
| 97 | |
---|
| 98 | def update(self): |
---|
| 99 | self._cache = {} |
---|
| 100 | |
---|
| 101 | def numpoints(self): |
---|
[7e224c2] | 102 | """ |
---|
| 103 | Return the number of points |
---|
| 104 | """ |
---|
[3c56da87] | 105 | return len(self.Iq) |
---|
[14de349] | 106 | |
---|
| 107 | def parameters(self): |
---|
[7e224c2] | 108 | """ |
---|
[346bc88] | 109 | Return a dictionary of parameters |
---|
[7e224c2] | 110 | """ |
---|
[346bc88] | 111 | return self.model.parameters() |
---|
[14de349] | 112 | |
---|
| 113 | def theory(self): |
---|
| 114 | if 'theory' not in self._cache: |
---|
[7cf2cfd] | 115 | pars = dict((k, v.value) for k,v in self.model.parameters().items()) |
---|
| 116 | self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff) |
---|
| 117 | """ |
---|
[abb22f4] | 118 | if self._fn is None: |
---|
[7cf2cfd] | 119 | q_input = self.model.kernel.make_input(self._kernel_inputs) |
---|
[346bc88] | 120 | self._fn = self.model.kernel(q_input) |
---|
[abb22f4] | 121 | |
---|
[346bc88] | 122 | fixed_pars = [getattr(self.model, p).value for p in self._fn.fixed_pars] |
---|
[abb22f4] | 123 | pd_pars = [self._get_weights(p) for p in self._fn.pd_pars] |
---|
[f1ecfa92] | 124 | #print fixed_pars,pd_pars |
---|
[346bc88] | 125 | Iq_calc = self._fn(fixed_pars, pd_pars, self.cutoff) |
---|
[abb22f4] | 126 | #self._theory[:] = self._fn.eval(pars, pd_pars) |
---|
[c97724e] | 127 | if self.data_type == 'sesans': |
---|
[3c56da87] | 128 | result = sesans.hankel(self.data.x, self.data.lam * 1e-9, |
---|
| 129 | self.data.sample.thickness / 10, |
---|
[7cf2cfd] | 130 | self._kernel_inputs[0], Iq_calc) |
---|
[3c56da87] | 131 | self._cache['theory'] = result |
---|
[c97724e] | 132 | else: |
---|
[346bc88] | 133 | Iq = self.resolution.apply(Iq_calc) |
---|
| 134 | self._cache['theory'] = Iq |
---|
[7cf2cfd] | 135 | """ |
---|
[14de349] | 136 | return self._cache['theory'] |
---|
| 137 | |
---|
| 138 | def residuals(self): |
---|
| 139 | #if np.any(self.err ==0): print "zeros in err" |
---|
[346bc88] | 140 | return (self.theory() - self.Iq) / self.dIq |
---|
[14de349] | 141 | |
---|
| 142 | def nllf(self): |
---|
[3c56da87] | 143 | delta = self.residuals() |
---|
[14de349] | 144 | #if np.any(np.isnan(R)): print "NaN in residuals" |
---|
[3c56da87] | 145 | return 0.5 * np.sum(delta ** 2) |
---|
[14de349] | 146 | |
---|
[3c56da87] | 147 | #def __call__(self): |
---|
| 148 | # return 2 * self.nllf() / self.dof |
---|
[14de349] | 149 | |
---|
| 150 | def plot(self, view='log'): |
---|
[c97724e] | 151 | """ |
---|
| 152 | Plot the data and residuals. |
---|
| 153 | """ |
---|
[7cf2cfd] | 154 | data, theory, resid = self._data, self.theory(), self.residuals() |
---|
| 155 | plot_theory(data, theory, resid, view) |
---|
[c97724e] | 156 | |
---|
| 157 | def simulate_data(self, noise=None): |
---|
[7cf2cfd] | 158 | Iq = self.theory() |
---|
| 159 | self._set_data(Iq, noise) |
---|
[14de349] | 160 | |
---|
| 161 | def save(self, basename): |
---|
| 162 | pass |
---|
| 163 | |
---|
[7cf2cfd] | 164 | def remove_get_weights(self, name): |
---|
[7e224c2] | 165 | """ |
---|
[aa4946b] | 166 | Get parameter dispersion weights |
---|
[7e224c2] | 167 | """ |
---|
[7cf2cfd] | 168 | info = self.model.kernel.info |
---|
| 169 | relative = name in info['partype']['pd-rel'] |
---|
| 170 | limits = info['limits'][name] |
---|
[3c56da87] | 171 | disperser, value, npts, width, nsigma = [ |
---|
[7cf2cfd] | 172 | getattr(self.model, name + ext) |
---|
[3c56da87] | 173 | for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')] |
---|
| 174 | value, weight = weights.get_weights( |
---|
[ce27e21] | 175 | disperser, int(npts.value), width.value, nsigma.value, |
---|
[7cf2cfd] | 176 | value.value, limits, relative) |
---|
[3c56da87] | 177 | return value, weight / np.sum(weight) |
---|
[ce27e21] | 178 | |
---|
[abb22f4] | 179 | def __getstate__(self): |
---|
| 180 | # Can't pickle gpu functions, so instead make them lazy |
---|
| 181 | state = self.__dict__.copy() |
---|
[7cf2cfd] | 182 | state['_kernel'] = None |
---|
[abb22f4] | 183 | return state |
---|
| 184 | |
---|
| 185 | def __setstate__(self, state): |
---|
[3c56da87] | 186 | # pylint: disable=attribute-defined-outside-init |
---|
[abb22f4] | 187 | self.__dict__ = state |
---|