source: sasmodels/sasmodels/weights.py @ d1fe925

gh-pages
Last change on this file since d1fe925 was d1fe925, checked in by ajj, 8 years ago

Creating gh_pages branch for docs

  • Property mode set to 100644
File size: 5.4 KB
Line 
1"""
2SAS distributions for polydispersity.
3"""
4# TODO: include dispersion docs with the disperser models
5from math import sqrt
6import numpy as np
7from scipy.special import gammaln
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        pars = {'type': self.type}
25        pars.update(self.__dict__)
26        return pars
27
28    # pylint: disable=no-self-use
29    def set_weights(self, values, weights):
30        raise RuntimeError("set_weights is only available for ArrayDispersion")
31
32    def get_weights(self, center, lb, ub, relative):
33        """
34        Return the weights for the distribution.
35
36        *center* is the center of the distribution
37
38        *lb*,*ub* are the min and max allowed values
39
40        *relative* is True if the distribution width is proportional to the
41        center value instead of absolute.  For polydispersity use relative.
42        For orientation parameters use absolute.
43        """
44        sigma = self.width * center if relative else self.width
45        if sigma == 0 or self.npts < 2:
46            if lb <= center <= ub:
47                return np.array([center], 'd'), np.array([1.], 'd')
48            else:
49                return np.array([], 'd'), np.array([], 'd')
50        return self._weights(center, sigma, lb, ub)
51
52    def _weights(self, center, sigma, lb, ub):
53        """actual work of computing the weights"""
54        raise NotImplementedError
55
56    def _linspace(self, center, sigma, lb, ub):
57        """helper function to provide linear spaced weight points within range"""
58        npts, nsigmas = self.npts, self.nsigmas
59        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts)
60        x = x[(x >= lb) & (x <= ub)]
61        return x
62
63
64class GaussianDispersion(Dispersion):
65    type = "gaussian"
66    default = dict(npts=35, width=0, nsigmas=3)
67    def _weights(self, center, sigma, lb, ub):
68        x = self._linspace(center, sigma, lb, ub)
69        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
70        return x, px
71
72
73class RectangleDispersion(Dispersion):
74    type = "rectangle"
75    default = dict(npts=35, width=0, nsigmas=1.70325)
76    def _weights(self, center, sigma, lb, ub):
77        x = self._linspace(center, sigma*sqrt(3.0), lb, ub)
78        px = np.ones_like(x)
79        return x, px
80
81
82class LogNormalDispersion(Dispersion):
83    type = "lognormal"
84    default = dict(npts=80, width=0, nsigmas=8)
85    def _weights(self, center, sigma, lb, ub):
86        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
87        px = np.exp(-0.5*(np.log(x)-center)**2)/sigma**2/(x*sigma)
88        return x, px
89
90
91class SchulzDispersion(Dispersion):
92    type = "schulz"
93    default = dict(npts=80, width=0, nsigmas=8)
94    def _weights(self, center, sigma, lb, ub):
95        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
96        R= x/center
97        z = (center/sigma)**2
98        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
99        px = np.exp(arg)
100        return x, px
101
102
103class ArrayDispersion(Dispersion):
104    type = "array"
105    default = dict(npts=35, width=0, nsigmas=1)
106    def __init__(self, npts=None, width=None, nsigmas=None):
107        Dispersion.__init__(self, npts, width, nsigmas)
108        self.values = np.array([0.], 'd')
109        self.weights = np.array([1.], 'd')
110
111    def set_weights(self, values, weights):
112        self.values = np.ascontiguousarray(values, 'd')
113        self.weights = np.ascontiguousarray(weights, 'd')
114        self.npts = len(values)
115
116    def _weights(self, center, sigma, lb, ub):
117        # TODO: interpolate into the array dispersion using npts
118        x = center + self.values*sigma
119        idx = (x>=lb)&(x<=ub)
120        x = x[idx]
121        px = self.weights[idx]
122        return x, px
123
124
125# dispersion name -> disperser lookup table.
126models = dict((d.type,d) for d in (
127    GaussianDispersion, RectangleDispersion,
128    ArrayDispersion, SchulzDispersion, LogNormalDispersion
129))
130
131
132def get_weights(disperser, n, width, nsigmas, value, limits, relative):
133    """
134    Return the set of values and weights for a polydisperse parameter.
135
136    *disperser* is the name of the disperser.
137
138    *n* is the number of points in the weight vector.
139
140    *width* is the width of the disperser distribution.
141
142    *nsigmas* is the number of sigmas to span for the dispersion convolution.
143
144    *value* is the value of the parameter in the model.
145
146    *limits* is [lb, ub], the lower and upper bound of the weight value.
147
148    *relative* is true if *width* is defined in proportion to the value
149    of the parameter, and false if it is an absolute width.
150
151    Returns *(value,weight)*, where *value* and *weight* are vectors.
152    """
153    cls = models[disperser]
154    obj = cls(n, width, nsigmas)
155    v,w = obj.get_weights(value, limits[0], limits[1], relative)
156    return v,w
157
158# Hack to allow sasview dispersion objects to interoperate with sasmodels
159dispersers = dict((v.__name__,k) for k,v in models.items())
160dispersers['DispersionModel'] = RectangleDispersion.type
161
Note: See TracBrowser for help on using the repository browser.