source: sasmodels/sasmodels/sasview_model.py @ d1fe925

gh-pages
Last change on this file since d1fe925 was d1fe925, checked in by ajj, 8 years ago

Creating gh_pages branch for docs

  • Property mode set to 100644
File size: 12.3 KB
Line 
1"""
2Sasview model constructor.
3
4Given a module defining an OpenCL kernel such as sasmodels.models.cylinder,
5create a sasview model class to run that kernel as follows::
6
7    from sasmodels.sasview_model import make_class
8    from sasmodels.models import cylinder
9    CylinderModel = make_class(cylinder, dtype='single')
10
11The model parameters for sasmodels are different from those in sasview.
12When reloading previously saved models, the parameters should be converted
13using :func:`sasmodels.convert.convert`.
14"""
15
16import math
17from copy import deepcopy
18import warnings
19
20import numpy as np
21
22from . import core
23
24def make_class(model_definition, dtype='single', namestyle='name'):
25    """
26    Load the sasview model defined in *kernel_module*.
27
28    Returns a class that can be used directly as a sasview model.
29
30    Defaults to using the new name for a model.  Setting
31    *namestyle='oldname'* will produce a class with a name
32    compatible with SasView.
33    """
34    model = core.load_model(model_definition, dtype=dtype)
35    def __init__(self, multfactor=1):
36        SasviewModel.__init__(self, model)
37    attrs = dict(__init__=__init__)
38    ConstructedModel = type(model.info[namestyle], (SasviewModel,), attrs)
39    return ConstructedModel
40
41class SasviewModel(object):
42    """
43    Sasview wrapper for opencl/ctypes model.
44    """
45    def __init__(self, model):
46        """Initialization"""
47        self._model = model
48
49        self.name = model.info['name']
50        self.oldname = model.info['oldname']
51        self.description = model.info['description']
52        self.category = None
53        self.multiplicity_info = None
54        self.is_multifunc = False
55
56        ## interpret the parameters
57        ## TODO: reorganize parameter handling
58        self.details = dict()
59        self.params = dict()
60        self.dispersion = dict()
61        partype = model.info['partype']
62        for name, units, default, limits, _, _ in model.info['parameters']:
63            self.params[name] = default
64            self.details[name] = [units] + limits
65
66        for name in partype['pd-2d']:
67            self.dispersion[name] = {
68                'width': 0,
69                'npts': 35,
70                'nsigmas': 3,
71                'type': 'gaussian',
72            }
73
74        self.orientation_params = (
75            partype['orientation']
76            + [n + '.width' for n in partype['orientation']]
77            + partype['magnetic'])
78        self.magnetic_params = partype['magnetic']
79        self.fixed = [n + '.width' for n in partype['pd-2d']]
80        self.non_fittable = []
81
82        ## independent parameter name and unit [string]
83        self.input_name = model.info.get("input_name", "Q")
84        self.input_unit = model.info.get("input_unit", "A^{-1}")
85        self.output_name = model.info.get("output_name", "Intensity")
86        self.output_unit = model.info.get("output_unit", "cm^{-1}")
87
88        ## _persistency_dict is used by sas.perspectives.fitting.basepage
89        ## to store dispersity reference.
90        ## TODO: _persistency_dict to persistency_dict throughout sasview
91        self._persistency_dict = {}
92
93        ## New fields introduced for opencl rewrite
94        self.cutoff = 1e-5
95
96    def __str__(self):
97        """
98        :return: string representation
99        """
100        return self.name
101
102    def is_fittable(self, par_name):
103        """
104        Check if a given parameter is fittable or not
105
106        :param par_name: the parameter name to check
107        """
108        return par_name.lower() in self.fixed
109        #For the future
110        #return self.params[str(par_name)].is_fittable()
111
112
113    # pylint: disable=no-self-use
114    def getProfile(self):
115        """
116        Get SLD profile
117
118        : return: (z, beta) where z is a list of depth of the transition points
119                beta is a list of the corresponding SLD values
120        """
121        return None, None
122
123    def setParam(self, name, value):
124        """
125        Set the value of a model parameter
126
127        :param name: name of the parameter
128        :param value: value of the parameter
129
130        """
131        # Look for dispersion parameters
132        toks = name.split('.')
133        if len(toks) == 2:
134            for item in self.dispersion.keys():
135                if item.lower() == toks[0].lower():
136                    for par in self.dispersion[item]:
137                        if par.lower() == toks[1].lower():
138                            self.dispersion[item][par] = value
139                            return
140        else:
141            # Look for standard parameter
142            for item in self.params.keys():
143                if item.lower() == name.lower():
144                    self.params[item] = value
145                    return
146
147        raise ValueError("Model does not contain parameter %s" % name)
148
149    def getParam(self, name):
150        """
151        Set the value of a model parameter
152
153        :param name: name of the parameter
154
155        """
156        # Look for dispersion parameters
157        toks = name.split('.')
158        if len(toks) == 2:
159            for item in self.dispersion.keys():
160                if item.lower() == toks[0].lower():
161                    for par in self.dispersion[item]:
162                        if par.lower() == toks[1].lower():
163                            return self.dispersion[item][par]
164        else:
165            # Look for standard parameter
166            for item in self.params.keys():
167                if item.lower() == name.lower():
168                    return self.params[item]
169
170        raise ValueError("Model does not contain parameter %s" % name)
171
172    def getParamList(self):
173        """
174        Return a list of all available parameters for the model
175        """
176        param_list = self.params.keys()
177        # WARNING: Extending the list with the dispersion parameters
178        param_list.extend(self.getDispParamList())
179        return param_list
180
181    def getDispParamList(self):
182        """
183        Return a list of all available parameters for the model
184        """
185        # TODO: fix test so that parameter order doesn't matter
186        ret = ['%s.%s' % (d.lower(), p)
187               for d in self._model.info['partype']['pd-2d']
188               for p in ('npts', 'nsigmas', 'width')]
189        #print ret
190        return ret
191
192    def clone(self):
193        """ Return a identical copy of self """
194        return deepcopy(self)
195
196    def run(self, x=0.0):
197        """
198        Evaluate the model
199
200        :param x: input q, or [q,phi]
201
202        :return: scattering function P(q)
203
204        **DEPRECATED**: use calculate_Iq instead
205        """
206        if isinstance(x, (list, tuple)):
207            # pylint: disable=unpacking-non-sequence
208            q, phi = x
209            return self.calculate_Iq([q * math.cos(phi)],
210                                     [q * math.sin(phi)])[0]
211        else:
212            return self.calculate_Iq([float(x)])[0]
213
214
215    def runXY(self, x=0.0):
216        """
217        Evaluate the model in cartesian coordinates
218
219        :param x: input q, or [qx, qy]
220
221        :return: scattering function P(q)
222
223        **DEPRECATED**: use calculate_Iq instead
224        """
225        if isinstance(x, (list, tuple)):
226            return self.calculate_Iq([float(x[0])], [float(x[1])])[0]
227        else:
228            return self.calculate_Iq([float(x)])[0]
229
230    def evalDistribution(self, qdist):
231        """
232        Evaluate a distribution of q-values.
233
234        * For 1D, a numpy array is expected as input: ::
235
236            evalDistribution(q)
237
238          where q is a numpy array.
239
240        * For 2D, a list of numpy arrays are expected: [qx,qy],
241          with 1D arrays::
242
243              qx = [ qx[0], qx[1], qx[2], ....]
244
245          and::
246
247              qy = [ qy[0], qy[1], qy[2], ....]
248
249        Then get ::
250
251            q = numpy.sqrt(qx^2+qy^2)
252
253        that is a qr in 1D array::
254
255            q = [q[0], q[1], q[2], ....]
256
257
258        :param qdist: ndarray of scalar q-values or list [qx,qy]
259        where qx,qy are 1D ndarrays
260        """
261        if isinstance(qdist, (list, tuple)):
262            # Check whether we have a list of ndarrays [qx,qy]
263            qx, qy = qdist
264            partype = self._model.info['partype']
265            if not partype['orientation'] and not partype['magnetic']:
266                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2))
267            else:
268                return self.calculate_Iq(qx, qy)
269
270        elif isinstance(qdist, np.ndarray):
271            # We have a simple 1D distribution of q-values
272            return self.calculate_Iq(qdist)
273
274        else:
275            raise TypeError("evalDistribution expects q or [qx, qy], not %r"
276                            % type(qdist))
277
278    def calculate_Iq(self, *args):
279        """
280        Calculate Iq for one set of q with the current parameters.
281
282        If the model is 1D, use *q*.  If 2D, use *qx*, *qy*.
283
284        This should NOT be used for fitting since it copies the *q* vectors
285        to the card for each evaluation.
286        """
287        q_vectors = [np.asarray(q) for q in args]
288        fn = self._model(self._model.make_input(q_vectors))
289        pars = [self.params[v] for v in fn.fixed_pars]
290        pd_pars = [self._get_weights(p) for p in fn.pd_pars]
291        result = fn(pars, pd_pars, self.cutoff)
292        fn.input.release()
293        fn.release()
294        return result
295
296    def calculate_ER(self):
297        """
298        Calculate the effective radius for P(q)*S(q)
299
300        :return: the value of the effective radius
301        """
302        ER = self._model.info.get('ER', None)
303        if ER is None:
304            return 1.0
305        else:
306            values, weights = self._dispersion_mesh()
307            fv = ER(*values)
308            #print values[0].shape, weights.shape, fv.shape
309            return np.sum(weights * fv) / np.sum(weights)
310
311    def calculate_VR(self):
312        """
313        Calculate the volf ratio for P(q)*S(q)
314
315        :return: the value of the volf ratio
316        """
317        VR = self._model.info.get('VR', None)
318        if VR is None:
319            return 1.0
320        else:
321            values, weights = self._dispersion_mesh()
322            whole, part = VR(*values)
323            return np.sum(weights * part) / np.sum(weights * whole)
324
325    def set_dispersion(self, parameter, dispersion):
326        """
327        Set the dispersion object for a model parameter
328
329        :param parameter: name of the parameter [string]
330        :param dispersion: dispersion object of type Dispersion
331        """
332        if parameter.lower() in (s.lower() for s in self.params.keys()):
333            # TODO: Store the disperser object directly in the model.
334            # The current method of creating one on the fly whenever it is
335            # needed is kind of funky.
336            # Note: can't seem to get disperser parameters from sasview
337            # (1) Could create a sasview model that has not yet # been
338            # converted, assign the disperser to one of its polydisperse
339            # parameters, then retrieve the disperser parameters from the
340            # sasview model.  (2) Could write a disperser parameter retriever
341            # in sasview.  (3) Could modify sasview to use sasmodels.weights
342            # dispersers.
343            # For now, rely on the fact that the sasview only ever uses
344            # new dispersers in the set_dispersion call and create a new
345            # one instead of trying to assign parameters.
346            from . import weights
347            disperser = weights.dispersers[dispersion.__class__.__name__]
348            dispersion = weights.models[disperser]()
349            self.dispersion[parameter] = dispersion.get_pars()
350        else:
351            raise ValueError("%r is not a dispersity or orientation parameter")
352
353    def _dispersion_mesh(self):
354        """
355        Create a mesh grid of dispersion parameters and weights.
356
357        Returns [p1,p2,...],w where pj is a vector of values for parameter j
358        and w is a vector containing the products for weights for each
359        parameter set in the vector.
360        """
361        pars = self._model.info['partype']['volume']
362        return core.dispersion_mesh([self._get_weights(p) for p in pars])
363
364    def _get_weights(self, par):
365        """
366            Return dispersion weights
367            :param par parameter name
368        """
369        from . import weights
370
371        relative = self._model.info['partype']['pd-rel']
372        limits = self._model.info['limits']
373        dis = self.dispersion[par]
374        value, weight = weights.get_weights(
375            dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
376            self.params[par], limits[par], par in relative)
377        return value, weight / np.sum(weight)
378
Note: See TracBrowser for help on using the repository browser.