source: sasmodels/sasmodels/weights.py

Last change on this file was d0b0f5d, checked in by GitHub <noreply@…>, 8 months ago

Merge pull request #67 from SasView?/ticket-608-user-defined-weights

Ticket 608 user defined weights

  • Property mode set to 100644
File size: 9.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
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
232SAS_WEIGHTS_PATH = "~/.sasview/weights"
233def load_weights(pattern=None):
234    # type: (str) -> None
235    """
236    Load dispersion distributions matching the given glob pattern
237    """
238    import logging
239    import os
240    import os.path
241    import glob
242    import traceback
243    from .custom import load_custom_kernel_module
244    if pattern is None:
245        path = os.environ.get("SAS_WEIGHTS_PATH", SAS_WEIGHTS_PATH)
246        pattern = os.path.join(path, "*.py")
247    for filename in sorted(glob.glob(os.path.expanduser(pattern))):
248        try:
249            #print("loading weights from", filename)
250            module = load_custom_kernel_module(filename)
251            MODELS[module.Dispersion.type] = module.Dispersion
252        except Exception as exc:
253            logging.error(traceback.format_exc(exc))
254
255def get_weights(disperser, n, width, nsigmas, value, limits, relative):
256    """
257    Return the set of values and weights for a polydisperse parameter.
258
259    *disperser* is the name of the disperser.
260
261    *n* is the number of points in the weight vector.
262
263    *width* is the width of the disperser distribution.
264
265    *nsigmas* is the number of sigmas to span for the dispersion convolution.
266
267    *value* is the value of the parameter in the model.
268
269    *limits* is [lb, ub], the lower and upper bound on the possible values.
270
271    *relative* is true if *width* is defined in proportion to the value
272    of the parameter, and false if it is an absolute width.
273
274    Returns *(value, weight)*, where *value* and *weight* are vectors.
275    """
276    if disperser == "array":
277        raise NotImplementedError("Don't handle arrays through get_weights;"
278                                  " use values and weights directly")
279    cls = MODELS[disperser]
280    obj = cls(n, width, nsigmas)
281    v, w = obj.get_weights(value, limits[0], limits[1], relative)
282    return v, w/np.sum(w)
283
284
285def plot_weights(model_info, mesh):
286    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None
287    """
288    Plot the weights returned by :func:`get_weights`.
289
290    *model_info* defines model parameters, etc.
291
292    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*)
293    for each parameter, where (*dispersity*, *weights*) pairs are the
294    distributions to be plotted.
295    """
296    import pylab
297
298    if any(len(dispersity) > 1 for value, dispersity, weights in mesh):
299        labels = [p.name for p in model_info.parameters.call_parameters]
300        #pylab.interactive(True)
301        pylab.figure()
302        for (v, x, w), s in zip(mesh, labels):
303            if len(x) > 1:
304                pylab.plot(x, w, '-o', label=s)
305        pylab.grid(True)
306        pylab.legend()
307        #pylab.show()
Note: See TracBrowser for help on using the repository browser.