source: sasmodels/sasmodels/weights.py @ e2592f0

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since e2592f0 was e2592f0, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

fix test warnings from py3 and numpy

  • Property mode set to 100644
File size: 9.0 KB
RevLine 
[ff7119b]1"""
2SAS distributions for polydispersity.
3"""
[ce27e21]4# TODO: include dispersion docs with the disperser models
[6cefbc9]5from __future__ import division, print_function
6
[7ae2b7f]7from math import sqrt  # type: ignore
[6cefbc9]8from collections import OrderedDict
9
[7ae2b7f]10import numpy as np  # type: ignore
11from scipy.special import gammaln  # type: ignore
[14de349]12
[733a3e1]13# TODO: include dispersion docs with the disperser models
14
[ce27e21]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):
[e2592f0]25        self.npts = self.default['npts'] if npts is None else int(npts)
[ce27e21]26        self.width = self.default['width'] if width is None else width
27        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas
[14de349]28
29    def get_pars(self):
[5c962df]30        """
31        Return the parameters to the disperser as a dictionary.
32        """
[ce27e21]33        pars = {'type': self.type}
34        pars.update(self.__dict__)
35        return pars
36
[3c56da87]37    # pylint: disable=no-self-use
[ce27e21]38    def set_weights(self, values, weights):
[5c962df]39        """
40        Set the weights on the disperser if it is :class:`ArrayDispersion`.
41        """
[ce27e21]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
[14de349]49
[823e620]50        *lb*, *ub* are the min and max allowed values
[ce27e21]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
[a5a12ca]57        if not relative:
58            # For orientation, the jitter is relative to 0 not the angle
59            center = 0
60            pass
[5d4777d]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')
[6cefbc9]66        x, px = self._weights(center, sigma, lb, ub)
67        return x, px
[ce27e21]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):
[5c962df]82    r"""
[0db85af]83    Gaussian dispersion, with 1-$\sigma$ width.
[5c962df]84
85    .. math::
86
87        w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right)
88    """
[ce27e21]89    type = "gaussian"
90    default = dict(npts=35, width=0, nsigmas=3)
91    def _weights(self, center, sigma, lb, ub):
[6cefbc9]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
[ce27e21]95        x = self._linspace(center, sigma, lb, ub)
[14de349]96        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
97        return x, px
98
[75e4319]99class UniformDispersion(Dispersion):
100    r"""
101    Uniform dispersion, with width $\sigma$.
102
103    .. math::
104
105        w = 1
106    """
107    type = "uniform"
[3d58247]108    default = dict(npts=35, width=0, nsigmas=None)
[75e4319]109    def _weights(self, center, sigma, lb, ub):
[3d58247]110        x = np.linspace(center-sigma, center+sigma, self.npts)
111        x = x[(x >= lb) & (x <= ub)]
[75e4319]112        return x, np.ones_like(x)
[ce27e21]113
114class RectangleDispersion(Dispersion):
[5c962df]115    r"""
116    Uniform dispersion, with width $\sqrt{3}\sigma$.
117
118    .. math::
119
120        w = 1
121    """
[ce27e21]122    type = "rectangle"
[1ceb951]123    default = dict(npts=35, width=0, nsigmas=1.73205)
[ce27e21]124    def _weights(self, center, sigma, lb, ub):
[3d58247]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)
[ce27e21]128
129class LogNormalDispersion(Dispersion):
[5c962df]130    r"""
[0db85af]131    log Gaussian dispersion, with 1-$\sigma$ width.
[5c962df]132
133    .. math::
134
135        w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma}
136    """
[ce27e21]137    type = "lognormal"
138    default = dict(npts=80, width=0, nsigmas=8)
139    def _weights(self, center, sigma, lb, ub):
[823e620]140        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
[f1a8811]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)
[ce27e21]144        return x, px
145
146
147class SchulzDispersion(Dispersion):
[5c962df]148    r"""
[0db85af]149    Schultz dispersion, with 1-$\sigma$ width.
[5c962df]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    """
[ce27e21]163    type = "schulz"
164    default = dict(npts=80, width=0, nsigmas=8)
165    def _weights(self, center, sigma, lb, ub):
[823e620]166        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8))
167        R = x/center
[ce27e21]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):
[5c962df]175    r"""
176    Empirical dispersion curve.
177
178    Use :meth:`set_weights` to set $w = f(x)$.
179    """
[ce27e21]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):
[5c962df]188        """
189        Set the weights for the given x values.
190        """
[ce27e21]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):
[6cefbc9]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
[823e620]200        idx = (x >= lb) & (x <= ub)
[ce27e21]201        x = x[idx]
202        px = self.weights[idx]
203        return x, px
204
[a5a12ca]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
[ce27e21]219
[ff7119b]220# dispersion name -> disperser lookup table.
[6cefbc9]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,
[75e4319]225    UniformDispersion,
[6cefbc9]226    ArrayDispersion,
227    LogNormalDispersion,
228    GaussianDispersion,
229    SchulzDispersion,
[a5a12ca]230    BoltzmannDispersion
[ce27e21]231))
232
[1780d59]233
234def get_weights(disperser, n, width, nsigmas, value, limits, relative):
[ff7119b]235    """
236    Return the set of values and weights for a polydisperse parameter.
237
238    *disperser* is the name of the disperser.
239
240    *n* is the number of points in the weight vector.
241
242    *width* is the width of the disperser distribution.
243
244    *nsigmas* is the number of sigmas to span for the dispersion convolution.
245
246    *value* is the value of the parameter in the model.
247
[6cefbc9]248    *limits* is [lb, ub], the lower and upper bound on the possible values.
[ff7119b]249
250    *relative* is true if *width* is defined in proportion to the value
251    of the parameter, and false if it is an absolute width.
252
[823e620]253    Returns *(value, weight)*, where *value* and *weight* are vectors.
[ff7119b]254    """
[6cefbc9]255    if disperser == "array":
[2d81cfe]256        raise NotImplementedError("Don't handle arrays through get_weights;"
257                                  " use values and weights directly")
[fa5fd8d]258    cls = MODELS[disperser]
[1780d59]259    obj = cls(n, width, nsigmas)
[823e620]260    v, w = obj.get_weights(value, limits[0], limits[1], relative)
[a5a12ca]261    return v, w/np.sum(w)
[6cefbc9]262
263
[a5a12ca]264def plot_weights(model_info, mesh):
265    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None
[6cefbc9]266    """
267    Plot the weights returned by :func:`get_weights`.
268
[a5a12ca]269    *model_info* defines model parameters, etc.
270
271    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*)
272    for each parameter, where (*dispersity*, *weights*) pairs are the
273    distributions to be plotted.
[6cefbc9]274    """
275    import pylab
276
[2d81cfe]277    if any(len(dispersity) > 1 for value, dispersity, weights in mesh):
[6cefbc9]278        labels = [p.name for p in model_info.parameters.call_parameters]
[a5a12ca]279        #pylab.interactive(True)
[6cefbc9]280        pylab.figure()
[2d81cfe]281        for (v, x, w), s in zip(mesh, labels):
[a5a12ca]282            if len(x) > 1:
283                pylab.plot(x, w, '-o', label=s)
[6cefbc9]284        pylab.grid(True)
285        pylab.legend()
[f1a8811]286        #pylab.show()
Note: See TracBrowser for help on using the repository browser.