source: sasmodels/sasmodels/weights.py @ 910c0f4

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

remove debugging trace

  • Property mode set to 100644
File size: 9.7 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
233SASMODELS_WEIGHTS = "~/.sasview/weights/*.py"
234def load_weights(pattern=None):
235    # type: (str) -> None
236    """
237    Load dispersion distributions matching the given 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        pattern = os.environ.get("SASMODELS_WEIGHTS", SASMODELS_WEIGHTS)
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.