source: sasmodels/sasmodels/weights.py @ 98c045a

ticket-1257-vesicle-product
Last change on this file since 98c045a was b297ba9, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

lint

  • Property mode set to 100644
File size: 9.0 KB
Line 
1"""
2SAS distributions for polydispersity.
3"""
4# TODO: include dispersion docs with the disperser models
5from __future__ import division, print_function
6
7from math import sqrt  # type: ignore
8from collections import OrderedDict
9
10import numpy as np  # type: ignore
11from scipy.special import gammaln  # type: ignore
12
13# TODO: include dispersion docs with the disperser models
14
15class Dispersion(object):
16    """
17    Base dispersion object.
18
19    Subclasses should define *_weights(center, sigma, lb, ub)*
20    which returns the x points and their corresponding weights.
21    """
22    type = "base disperser"
23    default = dict(npts=35, width=0, nsigmas=3)
24    def __init__(self, npts=None, width=None, nsigmas=None):
25        self.npts = self.default['npts'] if npts is None else int(npts)
26        self.width = self.default['width'] if width is None else width
27        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas
28
29    def get_pars(self):
30        """
31        Return the parameters to the disperser as a dictionary.
32        """
33        pars = {'type': self.type}
34        pars.update(self.__dict__)
35        return pars
36
37    # pylint: disable=no-self-use
38    def set_weights(self, values, weights):
39        """
40        Set the weights on the disperser if it is :class:`ArrayDispersion`.
41        """
42        raise RuntimeError("set_weights is only available for ArrayDispersion")
43
44    def get_weights(self, center, lb, ub, relative):
45        """
46        Return the weights for the distribution.
47
48        *center* is the center of the distribution
49
50        *lb*, *ub* are the min and max allowed values
51
52        *relative* is True if the distribution width is proportional to the
53        center value instead of absolute.  For polydispersity use relative.
54        For orientation parameters use absolute.
55        """
56        sigma = self.width * center if relative else self.width
57        if not relative:
58            # For orientation, the jitter is relative to 0 not the angle
59            center = 0
60        if sigma == 0 or self.npts < 2:
61            if lb <= center <= ub:
62                return np.array([center], 'd'), np.array([1.], 'd')
63            else:
64                return np.array([], 'd'), np.array([], 'd')
65        x, px = self._weights(center, sigma, lb, ub)
66        return x, px
67
68    def _weights(self, center, sigma, lb, ub):
69        """actual work of computing the weights"""
70        raise NotImplementedError
71
72    def _linspace(self, center, sigma, lb, ub):
73        """helper function to provide linear spaced weight points within range"""
74        npts, nsigmas = self.npts, self.nsigmas
75        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts)
76        x = x[(x >= lb) & (x <= ub)]
77        return x
78
79
80class GaussianDispersion(Dispersion):
81    r"""
82    Gaussian dispersion, with 1-$\sigma$ width.
83
84    .. math::
85
86        w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right)
87    """
88    type = "gaussian"
89    default = dict(npts=35, width=0, nsigmas=3)
90    def _weights(self, center, sigma, lb, ub):
91        # TODO: sample high probability regions more densely
92        # i.e., step uniformly in cumulative density rather than x value
93        # so weight = 1/Npts for all weights, but values are unevenly spaced
94        x = self._linspace(center, sigma, lb, ub)
95        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
96        return x, px
97
98class UniformDispersion(Dispersion):
99    r"""
100    Uniform dispersion, with width $\sigma$.
101
102    .. math::
103
104        w = 1
105    """
106    type = "uniform"
107    default = dict(npts=35, width=0, nsigmas=None)
108    def _weights(self, center, sigma, lb, ub):
109        x = np.linspace(center-sigma, center+sigma, self.npts)
110        x = x[(x >= lb) & (x <= ub)]
111        return x, np.ones_like(x)
112
113class RectangleDispersion(Dispersion):
114    r"""
115    Uniform dispersion, with width $\sqrt{3}\sigma$.
116
117    .. math::
118
119        w = 1
120    """
121    type = "rectangle"
122    default = dict(npts=35, width=0, nsigmas=1.73205)
123    def _weights(self, center, sigma, lb, ub):
124        x = self._linspace(center, sigma, lb, ub)
125        x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)]
126        return x, np.ones_like(x)
127
128class LogNormalDispersion(Dispersion):
129    r"""
130    log Gaussian dispersion, with 1-$\sigma$ width.
131
132    .. math::
133
134        w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma}
135    """
136    type = "lognormal"
137    default = dict(npts=80, width=0, nsigmas=8)
138    def _weights(self, center, sigma, lb, ub):
139        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
140        # sigma in the lognormal function is in ln(R) space, thus needs converting
141        sig = np.fabs(sigma/center)
142        px = np.exp(-0.5*((np.log(x)-np.log(center))/sig)**2)/(x*sig)
143        return x, px
144
145
146class SchulzDispersion(Dispersion):
147    r"""
148    Schultz dispersion, with 1-$\sigma$ width.
149
150    .. math::
151
152        w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)}
153
154    where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$.
155
156    This is evaluated using logarithms as
157
158    .. math::
159
160        w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right)
161    """
162    type = "schulz"
163    default = dict(npts=80, width=0, nsigmas=8)
164    def _weights(self, center, sigma, lb, ub):
165        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
166        R = x/center
167        z = (center/sigma)**2
168        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
169        px = np.exp(arg)
170        return x, px
171
172
173class ArrayDispersion(Dispersion):
174    r"""
175    Empirical dispersion curve.
176
177    Use :meth:`set_weights` to set $w = f(x)$.
178    """
179    type = "array"
180    default = dict(npts=35, width=0, nsigmas=1)
181    def __init__(self, npts=None, width=None, nsigmas=None):
182        Dispersion.__init__(self, npts, width, nsigmas)
183        self.values = np.array([0.], 'd')
184        self.weights = np.array([1.], 'd')
185
186    def set_weights(self, values, weights):
187        """
188        Set the weights for the given x values.
189        """
190        self.values = np.ascontiguousarray(values, 'd')
191        self.weights = np.ascontiguousarray(weights, 'd')
192        self.npts = len(values)
193
194    def _weights(self, center, sigma, lb, ub):
195        # TODO: rebin the array dispersion using npts
196        # TODO: use a distribution that can be recentered and scaled
197        x = self.values
198        #x = center + self.values*sigma
199        idx = (x >= lb) & (x <= ub)
200        x = x[idx]
201        px = self.weights[idx]
202        return x, px
203
204class BoltzmannDispersion(Dispersion):
205    r"""
206    Boltzmann dispersion, with $\sigma=k T/E$.
207
208    .. math::
209
210        w = \exp\left( -|x - c|/\sigma\right)
211    """
212    type = "boltzmann"
213    default = dict(npts=35, width=0, nsigmas=3)
214    def _weights(self, center, sigma, lb, ub):
215        x = self._linspace(center, sigma, lb, ub)
216        px = np.exp(-np.abs(x-center) / np.abs(sigma))
217        return x, px
218
219# dispersion name -> disperser lookup table.
220# Maintain order since this is used by sasview GUI to order the options in
221# the dispersion type combobox.
222MODELS = OrderedDict((d.type, d) for d in (
223    RectangleDispersion,
224    UniformDispersion,
225    ArrayDispersion,
226    LogNormalDispersion,
227    GaussianDispersion,
228    SchulzDispersion,
229    BoltzmannDispersion
230))
231
232
233def get_weights(disperser, n, width, nsigmas, value, limits, relative):
234    """
235    Return the set of values and weights for a polydisperse parameter.
236
237    *disperser* is the name of the disperser.
238
239    *n* is the number of points in the weight vector.
240
241    *width* is the width of the disperser distribution.
242
243    *nsigmas* is the number of sigmas to span for the dispersion convolution.
244
245    *value* is the value of the parameter in the model.
246
247    *limits* is [lb, ub], the lower and upper bound on the possible values.
248
249    *relative* is true if *width* is defined in proportion to the value
250    of the parameter, and false if it is an absolute width.
251
252    Returns *(value, weight)*, where *value* and *weight* are vectors.
253    """
254    if disperser == "array":
255        raise NotImplementedError("Don't handle arrays through get_weights;"
256                                  " use values and weights directly")
257    cls = MODELS[disperser]
258    obj = cls(n, width, nsigmas)
259    v, w = obj.get_weights(value, limits[0], limits[1], relative)
260    return v, w/np.sum(w)
261
262
263def plot_weights(model_info, mesh):
264    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None
265    """
266    Plot the weights returned by :func:`get_weights`.
267
268    *model_info* defines model parameters, etc.
269
270    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*)
271    for each parameter, where (*dispersity*, *weights*) pairs are the
272    distributions to be plotted.
273    """
274    import pylab
275
276    if any(len(dispersity) > 1 for value, dispersity, weights in mesh):
277        labels = [p.name for p in model_info.parameters.call_parameters]
278        #pylab.interactive(True)
279        pylab.figure()
280        for (v, x, w), s in zip(mesh, labels):
281            if len(x) > 1:
282                pylab.plot(x, w, '-o', label=s)
283        pylab.grid(True)
284        pylab.legend()
285        #pylab.show()
Note: See TracBrowser for help on using the repository browser.