source: sasmodels/sasmodels/weights.py @ 56b2687

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

Merge branch 'master' into polydisp

Conflicts:

README.rst
sasmodels/core.py
sasmodels/data.py
sasmodels/generate.py
sasmodels/kernelcl.py
sasmodels/kerneldll.py
sasmodels/sasview_model.py

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