source: sasmodels/sasmodels/weights.py @ 1780d59

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

hack sasview model polydispersity

  • Property mode set to 100644
File size: 4.5 KB
Line 
1# TODO: include dispersion docs with the disperser models
2from math import sqrt
3import numpy as np
4from scipy.special import gammaln
5
6class Dispersion(object):
7    """
8    Base dispersion object.
9
10    Subclasses should define *_weights(center, sigma, lb, ub)*
11    which returns the x points and their corresponding weights.
12    """
13    type = "base disperser"
14    default = dict(npts=35, width=0, nsigmas=3)
15    def __init__(self, npts=None, width=None, nsigmas=None):
16        self.npts = self.default['npts'] if npts is None else npts
17        self.width = self.default['width'] if width is None else width
18        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas
19
20    def get_pars(self):
21        pars = {'type': self.type}
22        pars.update(self.__dict__)
23        return pars
24
25    def set_weights(self, values, weights):
26        raise RuntimeError("set_weights is only available for ArrayDispersion")
27
28    def get_weights(self, center, lb, ub, relative):
29        """
30        Return the weights for the distribution.
31
32        *center* is the center of the distribution
33
34        *lb*,*ub* are the min and max allowed values
35
36        *relative* is True if the distribution width is proportional to the
37        center value instead of absolute.  For polydispersity use relative.
38        For orientation parameters use absolute.
39        """
40        sigma = self.width * center if relative else self.width
41        if sigma == 0:
42            return np.array([center], 'd'), np.array([1.], 'd')
43        return self._weights(center, sigma, lb, ub)
44
45    def _weights(self, center, sigma, lb, ub):
46        """actual work of computing the weights"""
47        raise NotImplementedError
48
49    def _linspace(self, center, sigma, lb, ub):
50        """helper function to provide linear spaced weight points within range"""
51        npts, nsigmas = self.npts, self.nsigmas
52        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts)
53        x = x[(x >= lb) & (x <= ub)]
54        return x
55
56
57class GaussianDispersion(Dispersion):
58    type = "gaussian"
59    default = dict(npts=35, width=0, nsigmas=3)
60    def _weights(self, center, sigma, lb, ub):
61        x = self._linspace(center, sigma, lb, ub)
62        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
63        return x, px
64
65
66class RectangleDispersion(Dispersion):
67    type = "rectangle"
68    default = dict(npts=35, width=0, nsigmas=1.70325)
69    def _weights(self, center, sigma, lb, ub):
70        x = self._linspace(center, sigma*sqrt(3.0), lb, ub)
71        px = np.ones_like(x)
72        return x, px
73
74
75class LogNormalDispersion(Dispersion):
76    type = "lognormal"
77    default = dict(npts=80, width=0, nsigmas=8)
78    def _weights(self, center, sigma, lb, ub):
79        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
80        px = np.exp(-0.5*(np.log(x)-center)**2)/sigma**2/(x*sigma)
81        return x, px
82
83
84class SchulzDispersion(Dispersion):
85    type = "schulz"
86    default = dict(npts=80, width=0, nsigmas=8)
87    def _weights(self, center, sigma, lb, ub):
88        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
89        R= x/center
90        z = (center/sigma)**2
91        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
92        px = np.exp(arg)
93        return x, px
94
95
96class ArrayDispersion(Dispersion):
97    type = "array"
98    default = dict(npts=35, width=0, nsigmas=1)
99    def __init__(self, npts=None, width=None, nsigmas=None):
100        Dispersion.__init__(self, npts, width, nsigmas)
101        self.values = np.array([0.], 'd')
102        self.weights = np.array([1.], 'd')
103
104    def set_weights(self, values, weights):
105        self.values = np.ascontiguousarray(values, 'd')
106        self.weights = np.ascontiguousarray(weights, 'd')
107        self.npts = len(values)
108
109    def _weights(self, center, sigma, lb, ub):
110        # TODO: interpolate into the array dispersion using npts
111        x = center + self.values*sigma
112        idx = (x>=lb)&(x<=ub)
113        x = x[idx]
114        px = self.weights[idx]
115        return x, px
116
117
118models = dict((d.type,d) for d in (
119    GaussianDispersion, RectangleDispersion,
120    ArrayDispersion, SchulzDispersion, LogNormalDispersion
121))
122
123
124def get_weights(disperser, n, width, nsigmas, value, limits, relative):
125    cls = models[disperser]
126    obj = cls(n, width, nsigmas)
127    v,w = obj.get_weights(value, limits[0], limits[1], relative)
128    return v,w
129
130# Hack to allow sasview dispersion objects to interoperate with sasmodels
131dispersers = dict((v.__name__,k) for k,v in models.items())
132dispersers['DispersionModel'] = RectangleDispersion.type
133
Note: See TracBrowser for help on using the repository browser.