[ff7119b] | 1 | """ |
---|
| 2 | SAS distributions for polydispersity. |
---|
| 3 | """ |
---|
[ce27e21] | 4 | # TODO: include dispersion docs with the disperser models |
---|
| 5 | from math import sqrt |
---|
[14de349] | 6 | import numpy as np |
---|
[ce27e21] | 7 | from scipy.special import gammaln |
---|
[14de349] | 8 | |
---|
[ce27e21] | 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 |
---|
[14de349] | 22 | |
---|
| 23 | def get_pars(self): |
---|
[5c962df] | 24 | """ |
---|
| 25 | Return the parameters to the disperser as a dictionary. |
---|
| 26 | """ |
---|
[ce27e21] | 27 | pars = {'type': self.type} |
---|
| 28 | pars.update(self.__dict__) |
---|
| 29 | return pars |
---|
| 30 | |
---|
[3c56da87] | 31 | # pylint: disable=no-self-use |
---|
[ce27e21] | 32 | def set_weights(self, values, weights): |
---|
[5c962df] | 33 | """ |
---|
| 34 | Set the weights on the disperser if it is :class:`ArrayDispersion`. |
---|
| 35 | """ |
---|
[ce27e21] | 36 | raise RuntimeError("set_weights is only available for ArrayDispersion") |
---|
| 37 | |
---|
| 38 | def get_weights(self, center, lb, ub, relative): |
---|
| 39 | """ |
---|
| 40 | Return the weights for the distribution. |
---|
| 41 | |
---|
| 42 | *center* is the center of the distribution |
---|
[14de349] | 43 | |
---|
[823e620] | 44 | *lb*, *ub* are the min and max allowed values |
---|
[ce27e21] | 45 | |
---|
| 46 | *relative* is True if the distribution width is proportional to the |
---|
| 47 | center value instead of absolute. For polydispersity use relative. |
---|
| 48 | For orientation parameters use absolute. |
---|
| 49 | """ |
---|
| 50 | sigma = self.width * center if relative else self.width |
---|
[5d4777d] | 51 | if sigma == 0 or self.npts < 2: |
---|
| 52 | if lb <= center <= ub: |
---|
| 53 | return np.array([center], 'd'), np.array([1.], 'd') |
---|
| 54 | else: |
---|
| 55 | return np.array([], 'd'), np.array([], 'd') |
---|
[ce27e21] | 56 | return self._weights(center, sigma, lb, ub) |
---|
| 57 | |
---|
| 58 | def _weights(self, center, sigma, lb, ub): |
---|
| 59 | """actual work of computing the weights""" |
---|
| 60 | raise NotImplementedError |
---|
| 61 | |
---|
| 62 | def _linspace(self, center, sigma, lb, ub): |
---|
| 63 | """helper function to provide linear spaced weight points within range""" |
---|
| 64 | npts, nsigmas = self.npts, self.nsigmas |
---|
| 65 | x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts) |
---|
| 66 | x = x[(x >= lb) & (x <= ub)] |
---|
| 67 | return x |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | class GaussianDispersion(Dispersion): |
---|
[5c962df] | 71 | r""" |
---|
| 72 | Gaussian dispersion, with 1-\ $\sigma$ width. |
---|
| 73 | |
---|
| 74 | .. math:: |
---|
| 75 | |
---|
| 76 | w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right) |
---|
| 77 | """ |
---|
[ce27e21] | 78 | type = "gaussian" |
---|
| 79 | default = dict(npts=35, width=0, nsigmas=3) |
---|
| 80 | def _weights(self, center, sigma, lb, ub): |
---|
| 81 | x = self._linspace(center, sigma, lb, ub) |
---|
[14de349] | 82 | px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) |
---|
| 83 | return x, px |
---|
| 84 | |
---|
[ce27e21] | 85 | |
---|
| 86 | class RectangleDispersion(Dispersion): |
---|
[5c962df] | 87 | r""" |
---|
| 88 | Uniform dispersion, with width $\sqrt{3}\sigma$. |
---|
| 89 | |
---|
| 90 | .. math:: |
---|
| 91 | |
---|
| 92 | w = 1 |
---|
| 93 | """ |
---|
[ce27e21] | 94 | type = "rectangle" |
---|
| 95 | default = dict(npts=35, width=0, nsigmas=1.70325) |
---|
| 96 | def _weights(self, center, sigma, lb, ub): |
---|
| 97 | x = self._linspace(center, sigma*sqrt(3.0), lb, ub) |
---|
| 98 | px = np.ones_like(x) |
---|
| 99 | return x, px |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | class LogNormalDispersion(Dispersion): |
---|
[5c962df] | 103 | r""" |
---|
| 104 | log Gaussian dispersion, with 1-\ $\sigma$ width. |
---|
| 105 | |
---|
| 106 | .. math:: |
---|
| 107 | |
---|
| 108 | w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma} |
---|
| 109 | """ |
---|
[ce27e21] | 110 | type = "lognormal" |
---|
| 111 | default = dict(npts=80, width=0, nsigmas=8) |
---|
| 112 | def _weights(self, center, sigma, lb, ub): |
---|
[823e620] | 113 | x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) |
---|
[5c962df] | 114 | px = np.exp(-0.5*(np.log(x)-center)**2/sigma**2)/(x*sigma) |
---|
[ce27e21] | 115 | return x, px |
---|
| 116 | |
---|
| 117 | |
---|
| 118 | class SchulzDispersion(Dispersion): |
---|
[5c962df] | 119 | r""" |
---|
| 120 | Schultz dispersion, with 1-\ $\sigma$ width. |
---|
| 121 | |
---|
| 122 | .. math:: |
---|
| 123 | |
---|
| 124 | w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)} |
---|
| 125 | |
---|
| 126 | where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$. |
---|
| 127 | |
---|
| 128 | This is evaluated using logarithms as |
---|
| 129 | |
---|
| 130 | .. math:: |
---|
| 131 | |
---|
| 132 | w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right) |
---|
| 133 | """ |
---|
[ce27e21] | 134 | type = "schulz" |
---|
| 135 | default = dict(npts=80, width=0, nsigmas=8) |
---|
| 136 | def _weights(self, center, sigma, lb, ub): |
---|
[823e620] | 137 | x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) |
---|
| 138 | R = x/center |
---|
[ce27e21] | 139 | z = (center/sigma)**2 |
---|
| 140 | arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z) |
---|
| 141 | px = np.exp(arg) |
---|
| 142 | return x, px |
---|
| 143 | |
---|
| 144 | |
---|
| 145 | class ArrayDispersion(Dispersion): |
---|
[5c962df] | 146 | r""" |
---|
| 147 | Empirical dispersion curve. |
---|
| 148 | |
---|
| 149 | Use :meth:`set_weights` to set $w = f(x)$. |
---|
| 150 | """ |
---|
[ce27e21] | 151 | type = "array" |
---|
| 152 | default = dict(npts=35, width=0, nsigmas=1) |
---|
| 153 | def __init__(self, npts=None, width=None, nsigmas=None): |
---|
| 154 | Dispersion.__init__(self, npts, width, nsigmas) |
---|
| 155 | self.values = np.array([0.], 'd') |
---|
| 156 | self.weights = np.array([1.], 'd') |
---|
| 157 | |
---|
| 158 | def set_weights(self, values, weights): |
---|
[5c962df] | 159 | """ |
---|
| 160 | Set the weights for the given x values. |
---|
| 161 | """ |
---|
[ce27e21] | 162 | self.values = np.ascontiguousarray(values, 'd') |
---|
| 163 | self.weights = np.ascontiguousarray(weights, 'd') |
---|
| 164 | self.npts = len(values) |
---|
| 165 | |
---|
| 166 | def _weights(self, center, sigma, lb, ub): |
---|
| 167 | # TODO: interpolate into the array dispersion using npts |
---|
| 168 | x = center + self.values*sigma |
---|
[823e620] | 169 | idx = (x >= lb) & (x <= ub) |
---|
[ce27e21] | 170 | x = x[idx] |
---|
| 171 | px = self.weights[idx] |
---|
| 172 | return x, px |
---|
| 173 | |
---|
| 174 | |
---|
[ff7119b] | 175 | # dispersion name -> disperser lookup table. |
---|
[823e620] | 176 | models = dict((d.type, d) for d in ( |
---|
[ce27e21] | 177 | GaussianDispersion, RectangleDispersion, |
---|
| 178 | ArrayDispersion, SchulzDispersion, LogNormalDispersion |
---|
| 179 | )) |
---|
| 180 | |
---|
[1780d59] | 181 | |
---|
| 182 | def get_weights(disperser, n, width, nsigmas, value, limits, relative): |
---|
[ff7119b] | 183 | """ |
---|
| 184 | Return the set of values and weights for a polydisperse parameter. |
---|
| 185 | |
---|
| 186 | *disperser* is the name of the disperser. |
---|
| 187 | |
---|
| 188 | *n* is the number of points in the weight vector. |
---|
| 189 | |
---|
| 190 | *width* is the width of the disperser distribution. |
---|
| 191 | |
---|
| 192 | *nsigmas* is the number of sigmas to span for the dispersion convolution. |
---|
| 193 | |
---|
| 194 | *value* is the value of the parameter in the model. |
---|
| 195 | |
---|
| 196 | *limits* is [lb, ub], the lower and upper bound of the weight value. |
---|
| 197 | |
---|
| 198 | *relative* is true if *width* is defined in proportion to the value |
---|
| 199 | of the parameter, and false if it is an absolute width. |
---|
| 200 | |
---|
[823e620] | 201 | Returns *(value, weight)*, where *value* and *weight* are vectors. |
---|
[ff7119b] | 202 | """ |
---|
[ce27e21] | 203 | cls = models[disperser] |
---|
[1780d59] | 204 | obj = cls(n, width, nsigmas) |
---|
[823e620] | 205 | v, w = obj.get_weights(value, limits[0], limits[1], relative) |
---|
| 206 | return v, w |
---|
[1780d59] | 207 | |
---|
| 208 | # Hack to allow sasview dispersion objects to interoperate with sasmodels |
---|
[823e620] | 209 | dispersers = dict((v.__name__, k) for k, v in models.items()) |
---|
[1780d59] | 210 | dispersers['DispersionModel'] = RectangleDispersion.type |
---|
| 211 | |
---|