source: sasmodels/sasmodels/weights.py @ f1a8811

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

conform to sasview definition of lognormal distribution

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