source: sasmodels/sasmodels/bumps_model.py @ cf404cb

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since cf404cb was 9404dd3, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

python 3.x support

  • Property mode set to 100644
File size: 6.2 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 the
9sasview data loader.  *Experiment* takes a *cutoff* parameter controlling
10how far the polydispersity integral extends.
11
12"""
13
14import datetime
15import warnings
16
17import numpy as np
18
19from . import sesans
20from . import weights
21from .data import plot_theory
22from .direct_model import DataMixin
23
24# CRUFT: old style bumps wrapper which doesn't separate data and model
25def 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
34class Model(object):
35    def __init__(self, model, **kw):
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
40        self._sasmodel = model
41        partype = model.info['partype']
42
43        pars = []
44        for p in model.info['parameters']:
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
73
74class Experiment(DataMixin):
75    """
76    Return a bumps wrapper for a SAS model.
77
78    *data* is the data to be fitted.
79
80    *model* is the SAS model from :func:`core.load_model`.
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    """
90    def __init__(self, data, model, cutoff=1e-5):
91
92        # remember inputs so we can inspect from outside
93        self.model = model
94        self.cutoff = cutoff
95        self._interpret_data(data, model._sasmodel)
96        self.update()
97
98    def update(self):
99        self._cache = {}
100
101    def numpoints(self):
102        """
103            Return the number of points
104        """
105        return len(self.Iq)
106
107    def parameters(self):
108        """
109        Return a dictionary of parameters
110        """
111        return self.model.parameters()
112
113    def theory(self):
114        if 'theory' not in self._cache:
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            """
118            if self._fn is None:
119                q_input = self.model.kernel.make_input(self._kernel_inputs)
120                self._fn = self.model.kernel(q_input)
121
122            fixed_pars = [getattr(self.model, p).value for p in self._fn.fixed_pars]
123            pd_pars = [self._get_weights(p) for p in self._fn.pd_pars]
124            #print(fixed_pars,pd_pars)
125            Iq_calc = self._fn(fixed_pars, pd_pars, self.cutoff)
126            #self._theory[:] = self._fn.eval(pars, pd_pars)
127            if self.data_type == 'sesans':
128                result = sesans.hankel(self.data.x, self.data.lam * 1e-9,
129                                       self.data.sample.thickness / 10,
130                                       self._kernel_inputs[0], Iq_calc)
131                self._cache['theory'] = result
132            else:
133                Iq = self.resolution.apply(Iq_calc)
134                self._cache['theory'] = Iq
135            """
136        return self._cache['theory']
137
138    def residuals(self):
139        #if np.any(self.err ==0): print("zeros in err")
140        return (self.theory() - self.Iq) / self.dIq
141
142    def nllf(self):
143        delta = self.residuals()
144        #if np.any(np.isnan(R)): print("NaN in residuals")
145        return 0.5 * np.sum(delta ** 2)
146
147    #def __call__(self):
148    #    return 2 * self.nllf() / self.dof
149
150    def plot(self, view='log'):
151        """
152        Plot the data and residuals.
153        """
154        data, theory, resid = self._data, self.theory(), self.residuals()
155        plot_theory(data, theory, resid, view)
156
157    def simulate_data(self, noise=None):
158        Iq = self.theory()
159        self._set_data(Iq, noise)
160
161    def save(self, basename):
162        pass
163
164    def remove_get_weights(self, name):
165        """
166        Get parameter dispersion weights
167        """
168        info = self.model.kernel.info
169        relative = name in info['partype']['pd-rel']
170        limits = info['limits'][name]
171        disperser, value, npts, width, nsigma = [
172            getattr(self.model, name + ext)
173            for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]
174        value, weight = weights.get_weights(
175            disperser, int(npts.value), width.value, nsigma.value,
176            value.value, limits, relative)
177        return value, weight / np.sum(weight)
178
179    def __getstate__(self):
180        # Can't pickle gpu functions, so instead make them lazy
181        state = self.__dict__.copy()
182        state['_kernel'] = None
183        return state
184
185    def __setstate__(self, state):
186        # pylint: disable=attribute-defined-outside-init
187        self.__dict__ = state
Note: See TracBrowser for help on using the repository browser.