source: sasmodels/sasmodels/weights.py @ 733a3e1

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

minor code cleanup

  • Property mode set to 100644
File size: 6.4 KB
Line 
1"""
2SAS distributions for polydispersity.
3"""
4from __future__ import division
5from math import sqrt  # type: ignore
6import numpy as np  # type: ignore
7from scipy.special import gammaln  # type: ignore
8
9# TODO: include dispersion docs with the disperser models
10
11class Dispersion(object):
12    """
13    Base dispersion object.
14
15    Subclasses should define *_weights(center, sigma, lb, ub)*
16    which returns the x points and their corresponding weights.
17    """
18    type = "base disperser"
19    default = dict(npts=35, width=0, nsigmas=3)
20    def __init__(self, npts=None, width=None, nsigmas=None):
21        self.npts = self.default['npts'] if npts is None else npts
22        self.width = self.default['width'] if width is None else width
23        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas
24
25    def get_pars(self):
26        """
27        Return the parameters to the disperser as a dictionary.
28        """
29        pars = {'type': self.type}
30        pars.update(self.__dict__)
31        return pars
32
33    # pylint: disable=no-self-use
34    def set_weights(self, values, weights):
35        """
36        Set the weights on the disperser if it is :class:`ArrayDispersion`.
37        """
38        raise RuntimeError("set_weights is only available for ArrayDispersion")
39
40    def get_weights(self, center, lb, ub, relative):
41        """
42        Return the weights for the distribution.
43
44        *center* is the center of the distribution
45
46        *lb*, *ub* are the min and max allowed values
47
48        *relative* is True if the distribution width is proportional to the
49        center value instead of absolute.  For polydispersity use relative.
50        For orientation parameters use absolute.
51        """
52        sigma = self.width * center if relative else self.width
53        if sigma == 0 or self.npts < 2:
54            if lb <= center <= ub:
55                return np.array([center], 'd'), np.array([1.], 'd')
56            else:
57                return np.array([], 'd'), np.array([], 'd')
58        return self._weights(center, sigma, lb, ub)
59
60    def _weights(self, center, sigma, lb, ub):
61        """actual work of computing the weights"""
62        raise NotImplementedError
63
64    def _linspace(self, center, sigma, lb, ub):
65        """helper function to provide linear spaced weight points within range"""
66        npts, nsigmas = self.npts, self.nsigmas
67        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts)
68        x = x[(x >= lb) & (x <= ub)]
69        return x
70
71
72class GaussianDispersion(Dispersion):
73    r"""
74    Gaussian dispersion, with 1-\ $\sigma$ width.
75
76    .. math::
77
78        w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right)
79    """
80    type = "gaussian"
81    default = dict(npts=35, width=0, nsigmas=3)
82    def _weights(self, center, sigma, lb, ub):
83        x = self._linspace(center, sigma, lb, ub)
84        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
85        return x, px
86
87
88class RectangleDispersion(Dispersion):
89    r"""
90    Uniform dispersion, with width $\sqrt{3}\sigma$.
91
92    .. math::
93
94        w = 1
95    """
96    type = "rectangle"
97    default = dict(npts=35, width=0, nsigmas=1.70325)
98    def _weights(self, center, sigma, lb, ub):
99        x = self._linspace(center, sigma*sqrt(3.0), lb, ub)
100        px = np.ones_like(x)
101        return x, px
102
103
104class LogNormalDispersion(Dispersion):
105    r"""
106    log Gaussian dispersion, with 1-\ $\sigma$ width.
107
108    .. math::
109
110        w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma}
111    """
112    type = "lognormal"
113    default = dict(npts=80, width=0, nsigmas=8)
114    def _weights(self, center, sigma, lb, ub):
115        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
116        px = np.exp(-0.5*(np.log(x)-center)**2/sigma**2)/(x*sigma)
117        return x, px
118
119
120class SchulzDispersion(Dispersion):
121    r"""
122    Schultz dispersion, with 1-\ $\sigma$ width.
123
124    .. math::
125
126        w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)}
127
128    where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$.
129
130    This is evaluated using logarithms as
131
132    .. math::
133
134        w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right)
135    """
136    type = "schulz"
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        R = x/center
141        z = (center/sigma)**2
142        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
143        px = np.exp(arg)
144        return x, px
145
146
147class ArrayDispersion(Dispersion):
148    r"""
149    Empirical dispersion curve.
150
151    Use :meth:`set_weights` to set $w = f(x)$.
152    """
153    type = "array"
154    default = dict(npts=35, width=0, nsigmas=1)
155    def __init__(self, npts=None, width=None, nsigmas=None):
156        Dispersion.__init__(self, npts, width, nsigmas)
157        self.values = np.array([0.], 'd')
158        self.weights = np.array([1.], 'd')
159
160    def set_weights(self, values, weights):
161        """
162        Set the weights for the given x values.
163        """
164        self.values = np.ascontiguousarray(values, 'd')
165        self.weights = np.ascontiguousarray(weights, 'd')
166        self.npts = len(values)
167
168    def _weights(self, center, sigma, lb, ub):
169        # TODO: interpolate into the array dispersion using npts
170        x = center + self.values*sigma
171        idx = (x >= lb) & (x <= ub)
172        x = x[idx]
173        px = self.weights[idx]
174        return x, px
175
176
177# dispersion name -> disperser lookup table.
178MODELS = dict((d.type, d) for d in (
179    GaussianDispersion, RectangleDispersion,
180    ArrayDispersion, SchulzDispersion, LogNormalDispersion
181))
182
183
184def get_weights(disperser, n, width, nsigmas, value, limits, relative):
185    """
186    Return the set of values and weights for a polydisperse parameter.
187
188    *disperser* is the name of the disperser.
189
190    *n* is the number of points in the weight vector.
191
192    *width* is the width of the disperser distribution.
193
194    *nsigmas* is the number of sigmas to span for the dispersion convolution.
195
196    *value* is the value of the parameter in the model.
197
198    *limits* is [lb, ub], the lower and upper bound of the weight value.
199
200    *relative* is true if *width* is defined in proportion to the value
201    of the parameter, and false if it is an absolute width.
202
203    Returns *(value, weight)*, where *value* and *weight* are vectors.
204    """
205    cls = MODELS[disperser]
206    obj = cls(n, width, nsigmas)
207    v, w = obj.get_weights(value, limits[0], limits[1], relative)
208    return v, w
Note: See TracBrowser for help on using the repository browser.