source: sasmodels/sasmodels/weights.py @ f41027b

Last change on this file since f41027b was f41027b, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

change environment variable for distributions to SAS_WEIGHTS_PATH

  • 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 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            pass
61        if sigma == 0 or self.npts < 2:
62            if lb <= center <= ub:
63                return np.array([center], 'd'), np.array([1.], 'd')
64            else:
65                return np.array([], 'd'), np.array([], 'd')
66        x, px = self._weights(center, sigma, lb, ub)
67        return x, px
68
69    def _weights(self, center, sigma, lb, ub):
70        """actual work of computing the weights"""
71        raise NotImplementedError
72
73    def _linspace(self, center, sigma, lb, ub):
74        """helper function to provide linear spaced weight points within range"""
75        npts, nsigmas = self.npts, self.nsigmas
76        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts)
77        x = x[(x >= lb) & (x <= ub)]
78        return x
79
80
81class GaussianDispersion(Dispersion):
82    r"""
83    Gaussian dispersion, with 1-$\sigma$ width.
84
85    .. math::
86
87        w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right)
88    """
89    type = "gaussian"
90    default = dict(npts=35, width=0, nsigmas=3)
91    def _weights(self, center, sigma, lb, ub):
92        # TODO: sample high probability regions more densely
93        # i.e., step uniformly in cumulative density rather than x value
94        # so weight = 1/Npts for all weights, but values are unevenly spaced
95        x = self._linspace(center, sigma, lb, ub)
96        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
97        return x, px
98
99class UniformDispersion(Dispersion):
100    r"""
101    Uniform dispersion, with width $\sigma$.
102
103    .. math::
104
105        w = 1
106    """
107    type = "uniform"
108    default = dict(npts=35, width=0, nsigmas=None)
109    def _weights(self, center, sigma, lb, ub):
110        x = np.linspace(center-sigma, center+sigma, self.npts)
111        x = x[(x >= lb) & (x <= ub)]
112        return x, np.ones_like(x)
113
114class RectangleDispersion(Dispersion):
115    r"""
116    Uniform dispersion, with width $\sqrt{3}\sigma$.
117
118    .. math::
119
120        w = 1
121    """
122    type = "rectangle"
123    default = dict(npts=35, width=0, nsigmas=1.73205)
124    def _weights(self, center, sigma, lb, ub):
125         x = self._linspace(center, sigma, lb, ub)
126         x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)]
127         return x, np.ones_like(x)
128
129class LogNormalDispersion(Dispersion):
130    r"""
131    log Gaussian dispersion, with 1-$\sigma$ width.
132
133    .. math::
134
135        w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma}
136    """
137    type = "lognormal"
138    default = dict(npts=80, width=0, nsigmas=8)
139    def _weights(self, center, sigma, lb, ub):
140        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
141        # sigma in the lognormal function is in ln(R) space, thus needs converting
142        sig = np.fabs(sigma/center)
143        px = np.exp(-0.5*((np.log(x)-np.log(center))/sig)**2)/(x*sig)
144        return x, px
145
146
147class SchulzDispersion(Dispersion):
148    r"""
149    Schultz dispersion, with 1-$\sigma$ width.
150
151    .. math::
152
153        w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)}
154
155    where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$.
156
157    This is evaluated using logarithms as
158
159    .. math::
160
161        w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right)
162    """
163    type = "schulz"
164    default = dict(npts=80, width=0, nsigmas=8)
165    def _weights(self, center, sigma, lb, ub):
166        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
167        R = x/center
168        z = (center/sigma)**2
169        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
170        px = np.exp(arg)
171        return x, px
172
173
174class ArrayDispersion(Dispersion):
175    r"""
176    Empirical dispersion curve.
177
178    Use :meth:`set_weights` to set $w = f(x)$.
179    """
180    type = "array"
181    default = dict(npts=35, width=0, nsigmas=1)
182    def __init__(self, npts=None, width=None, nsigmas=None):
183        Dispersion.__init__(self, npts, width, nsigmas)
184        self.values = np.array([0.], 'd')
185        self.weights = np.array([1.], 'd')
186
187    def set_weights(self, values, weights):
188        """
189        Set the weights for the given x values.
190        """
191        self.values = np.ascontiguousarray(values, 'd')
192        self.weights = np.ascontiguousarray(weights, 'd')
193        self.npts = len(values)
194
195    def _weights(self, center, sigma, lb, ub):
196        # TODO: rebin the array dispersion using npts
197        # TODO: use a distribution that can be recentered and scaled
198        x = self.values
199        #x = center + self.values*sigma
200        idx = (x >= lb) & (x <= ub)
201        x = x[idx]
202        px = self.weights[idx]
203        return x, px
204
205class BoltzmannDispersion(Dispersion):
206    r"""
207    Boltzmann dispersion, with $\sigma=k T/E$.
208
209    .. math::
210
211        w = \exp\left( -|x - c|/\sigma\right)
212    """
213    type = "boltzmann"
214    default = dict(npts=35, width=0, nsigmas=3)
215    def _weights(self, center, sigma, lb, ub):
216        x = self._linspace(center, sigma, lb, ub)
217        px = np.exp(-np.abs(x-center) / np.abs(sigma))
218        return x, px
219
220# dispersion name -> disperser lookup table.
221# Maintain order since this is used by sasview GUI to order the options in
222# the dispersion type combobox.
223MODELS = OrderedDict((d.type, d) for d in (
224    RectangleDispersion,
225    UniformDispersion,
226    ArrayDispersion,
227    LogNormalDispersion,
228    GaussianDispersion,
229    SchulzDispersion,
230    BoltzmannDispersion
231))
232
233SAS_WEIGHTS_PATH = "~/.sasview/weights"
234def load_weights(pattern=None):
235    # type: (str) -> None
236    """
237    Load dispersion distributions matching the given glob pattern
238    """
239    import logging
240    import os
241    import os.path
242    import glob
243    import traceback
244    from .custom import load_custom_kernel_module
245    if pattern is None:
246        path = os.environ.get("SAS_WEIGHTS_PATH", SAS_WEIGHTS_PATH)
247        pattern = os.path.join(path, "*.py")
248    for filename in sorted(glob.glob(os.path.expanduser(pattern))):
249        try:
250            #print("loading weights from", filename)
251            module = load_custom_kernel_module(filename)
252            MODELS[module.Dispersion.type] = module.Dispersion
253        except Exception as exc:
254            logging.error(traceback.format_exc(exc))
255
256def get_weights(disperser, n, width, nsigmas, value, limits, relative):
257    """
258    Return the set of values and weights for a polydisperse parameter.
259
260    *disperser* is the name of the disperser.
261
262    *n* is the number of points in the weight vector.
263
264    *width* is the width of the disperser distribution.
265
266    *nsigmas* is the number of sigmas to span for the dispersion convolution.
267
268    *value* is the value of the parameter in the model.
269
270    *limits* is [lb, ub], the lower and upper bound on the possible values.
271
272    *relative* is true if *width* is defined in proportion to the value
273    of the parameter, and false if it is an absolute width.
274
275    Returns *(value, weight)*, where *value* and *weight* are vectors.
276    """
277    if disperser == "array":
278        raise NotImplementedError("Don't handle arrays through get_weights;"
279                                  " use values and weights directly")
280    cls = MODELS[disperser]
281    obj = cls(n, width, nsigmas)
282    v, w = obj.get_weights(value, limits[0], limits[1], relative)
283    return v, w/np.sum(w)
284
285
286def plot_weights(model_info, mesh):
287    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None
288    """
289    Plot the weights returned by :func:`get_weights`.
290
291    *model_info* defines model parameters, etc.
292
293    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*)
294    for each parameter, where (*dispersity*, *weights*) pairs are the
295    distributions to be plotted.
296    """
297    import pylab
298
299    if any(len(dispersity) > 1 for value, dispersity, weights in mesh):
300        labels = [p.name for p in model_info.parameters.call_parameters]
301        #pylab.interactive(True)
302        pylab.figure()
303        for (v, x, w), s in zip(mesh, labels):
304            if len(x) > 1:
305                pylab.plot(x, w, '-o', label=s)
306        pylab.grid(True)
307        pylab.legend()
308        #pylab.show()
Note: See TracBrowser for help on using the repository browser.