source: sasmodels/sasmodels/sasview_model.py @ 87c722e

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 87c722e was 87c722e, checked in by pkienzle, 9 years ago

use sas instead of sans

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