1 | """ |
---|
2 | SAS distributions for polydispersity. |
---|
3 | """ |
---|
4 | # TODO: include dispersion docs with the disperser models |
---|
5 | from math import sqrt |
---|
6 | import numpy as np |
---|
7 | from scipy.special import gammaln |
---|
8 | |
---|
9 | class 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 | |
---|
64 | class 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 | |
---|
73 | class 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 | |
---|
82 | class 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 | |
---|
91 | class 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 | |
---|
103 | class 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. |
---|
126 | models = dict((d.type,d) for d in ( |
---|
127 | GaussianDispersion, RectangleDispersion, |
---|
128 | ArrayDispersion, SchulzDispersion, LogNormalDispersion |
---|
129 | )) |
---|
130 | |
---|
131 | |
---|
132 | def 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 |
---|
159 | dispersers = dict((v.__name__,k) for k,v in models.items()) |
---|
160 | dispersers['DispersionModel'] = RectangleDispersion.type |
---|
161 | |
---|