[ff7119b] | 1 | """ |
---|
| 2 | SAS distributions for polydispersity. |
---|
| 3 | """ |
---|
[ce27e21] | 4 | # TODO: include dispersion docs with the disperser models |
---|
[6cefbc9] | 5 | from __future__ import division, print_function |
---|
| 6 | |
---|
[7ae2b7f] | 7 | from math import sqrt # type: ignore |
---|
[6cefbc9] | 8 | from collections import OrderedDict |
---|
| 9 | |
---|
[7ae2b7f] | 10 | import numpy as np # type: ignore |
---|
| 11 | from scipy.special import gammaln # type: ignore |
---|
[14de349] | 12 | |
---|
[733a3e1] | 13 | # TODO: include dispersion docs with the disperser models |
---|
| 14 | |
---|
[ce27e21] | 15 | class 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 |
---|
[5d4777d] | 60 | if sigma == 0 or self.npts < 2: |
---|
| 61 | if lb <= center <= ub: |
---|
| 62 | return np.array([center], 'd'), np.array([1.], 'd') |
---|
| 63 | else: |
---|
| 64 | return np.array([], 'd'), np.array([], 'd') |
---|
[6cefbc9] | 65 | x, px = self._weights(center, sigma, lb, ub) |
---|
| 66 | return x, px |
---|
[ce27e21] | 67 | |
---|
| 68 | def _weights(self, center, sigma, lb, ub): |
---|
| 69 | """actual work of computing the weights""" |
---|
| 70 | raise NotImplementedError |
---|
| 71 | |
---|
| 72 | def _linspace(self, center, sigma, lb, ub): |
---|
| 73 | """helper function to provide linear spaced weight points within range""" |
---|
| 74 | npts, nsigmas = self.npts, self.nsigmas |
---|
| 75 | x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts) |
---|
| 76 | x = x[(x >= lb) & (x <= ub)] |
---|
| 77 | return x |
---|
| 78 | |
---|
| 79 | |
---|
| 80 | class GaussianDispersion(Dispersion): |
---|
[5c962df] | 81 | r""" |
---|
[0db85af] | 82 | Gaussian dispersion, with 1-$\sigma$ width. |
---|
[5c962df] | 83 | |
---|
| 84 | .. math:: |
---|
| 85 | |
---|
| 86 | w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right) |
---|
| 87 | """ |
---|
[ce27e21] | 88 | type = "gaussian" |
---|
| 89 | default = dict(npts=35, width=0, nsigmas=3) |
---|
| 90 | def _weights(self, center, sigma, lb, ub): |
---|
[6cefbc9] | 91 | # TODO: sample high probability regions more densely |
---|
| 92 | # i.e., step uniformly in cumulative density rather than x value |
---|
| 93 | # so weight = 1/Npts for all weights, but values are unevenly spaced |
---|
[ce27e21] | 94 | x = self._linspace(center, sigma, lb, ub) |
---|
[14de349] | 95 | px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) |
---|
| 96 | return x, px |
---|
| 97 | |
---|
[75e4319] | 98 | class UniformDispersion(Dispersion): |
---|
| 99 | r""" |
---|
| 100 | Uniform dispersion, with width $\sigma$. |
---|
| 101 | |
---|
| 102 | .. math:: |
---|
| 103 | |
---|
| 104 | w = 1 |
---|
| 105 | """ |
---|
| 106 | type = "uniform" |
---|
[3d58247] | 107 | default = dict(npts=35, width=0, nsigmas=None) |
---|
[75e4319] | 108 | def _weights(self, center, sigma, lb, ub): |
---|
[3d58247] | 109 | x = np.linspace(center-sigma, center+sigma, self.npts) |
---|
| 110 | x = x[(x >= lb) & (x <= ub)] |
---|
[75e4319] | 111 | return x, np.ones_like(x) |
---|
[ce27e21] | 112 | |
---|
| 113 | class RectangleDispersion(Dispersion): |
---|
[5c962df] | 114 | r""" |
---|
| 115 | Uniform dispersion, with width $\sqrt{3}\sigma$. |
---|
| 116 | |
---|
| 117 | .. math:: |
---|
| 118 | |
---|
| 119 | w = 1 |
---|
| 120 | """ |
---|
[ce27e21] | 121 | type = "rectangle" |
---|
[1ceb951] | 122 | default = dict(npts=35, width=0, nsigmas=1.73205) |
---|
[ce27e21] | 123 | def _weights(self, center, sigma, lb, ub): |
---|
[b297ba9] | 124 | x = self._linspace(center, sigma, lb, ub) |
---|
| 125 | x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] |
---|
| 126 | return x, np.ones_like(x) |
---|
[ce27e21] | 127 | |
---|
| 128 | class LogNormalDispersion(Dispersion): |
---|
[5c962df] | 129 | r""" |
---|
[0db85af] | 130 | log Gaussian dispersion, with 1-$\sigma$ width. |
---|
[5c962df] | 131 | |
---|
| 132 | .. math:: |
---|
| 133 | |
---|
| 134 | w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma} |
---|
| 135 | """ |
---|
[ce27e21] | 136 | type = "lognormal" |
---|
| 137 | default = dict(npts=80, width=0, nsigmas=8) |
---|
| 138 | def _weights(self, center, sigma, lb, ub): |
---|
[823e620] | 139 | x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) |
---|
[f1a8811] | 140 | # sigma in the lognormal function is in ln(R) space, thus needs converting |
---|
| 141 | sig = np.fabs(sigma/center) |
---|
| 142 | px = np.exp(-0.5*((np.log(x)-np.log(center))/sig)**2)/(x*sig) |
---|
[ce27e21] | 143 | return x, px |
---|
| 144 | |
---|
| 145 | |
---|
| 146 | class SchulzDispersion(Dispersion): |
---|
[5c962df] | 147 | r""" |
---|
[0db85af] | 148 | Schultz dispersion, with 1-$\sigma$ width. |
---|
[5c962df] | 149 | |
---|
| 150 | .. math:: |
---|
| 151 | |
---|
| 152 | w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)} |
---|
| 153 | |
---|
| 154 | where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$. |
---|
| 155 | |
---|
| 156 | This is evaluated using logarithms as |
---|
| 157 | |
---|
| 158 | .. math:: |
---|
| 159 | |
---|
| 160 | w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right) |
---|
| 161 | """ |
---|
[ce27e21] | 162 | type = "schulz" |
---|
| 163 | default = dict(npts=80, width=0, nsigmas=8) |
---|
| 164 | def _weights(self, center, sigma, lb, ub): |
---|
[823e620] | 165 | x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) |
---|
| 166 | R = x/center |
---|
[ce27e21] | 167 | z = (center/sigma)**2 |
---|
| 168 | arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z) |
---|
| 169 | px = np.exp(arg) |
---|
| 170 | return x, px |
---|
| 171 | |
---|
| 172 | |
---|
| 173 | class ArrayDispersion(Dispersion): |
---|
[5c962df] | 174 | r""" |
---|
| 175 | Empirical dispersion curve. |
---|
| 176 | |
---|
| 177 | Use :meth:`set_weights` to set $w = f(x)$. |
---|
| 178 | """ |
---|
[ce27e21] | 179 | type = "array" |
---|
| 180 | default = dict(npts=35, width=0, nsigmas=1) |
---|
| 181 | def __init__(self, npts=None, width=None, nsigmas=None): |
---|
| 182 | Dispersion.__init__(self, npts, width, nsigmas) |
---|
| 183 | self.values = np.array([0.], 'd') |
---|
| 184 | self.weights = np.array([1.], 'd') |
---|
| 185 | |
---|
| 186 | def set_weights(self, values, weights): |
---|
[5c962df] | 187 | """ |
---|
| 188 | Set the weights for the given x values. |
---|
| 189 | """ |
---|
[ce27e21] | 190 | self.values = np.ascontiguousarray(values, 'd') |
---|
| 191 | self.weights = np.ascontiguousarray(weights, 'd') |
---|
| 192 | self.npts = len(values) |
---|
| 193 | |
---|
| 194 | def _weights(self, center, sigma, lb, ub): |
---|
[6cefbc9] | 195 | # TODO: rebin the array dispersion using npts |
---|
| 196 | # TODO: use a distribution that can be recentered and scaled |
---|
| 197 | x = self.values |
---|
| 198 | #x = center + self.values*sigma |
---|
[823e620] | 199 | idx = (x >= lb) & (x <= ub) |
---|
[ce27e21] | 200 | x = x[idx] |
---|
| 201 | px = self.weights[idx] |
---|
| 202 | return x, px |
---|
| 203 | |
---|
[a5a12ca] | 204 | class BoltzmannDispersion(Dispersion): |
---|
| 205 | r""" |
---|
| 206 | Boltzmann dispersion, with $\sigma=k T/E$. |
---|
| 207 | |
---|
| 208 | .. math:: |
---|
| 209 | |
---|
| 210 | w = \exp\left( -|x - c|/\sigma\right) |
---|
| 211 | """ |
---|
| 212 | type = "boltzmann" |
---|
| 213 | default = dict(npts=35, width=0, nsigmas=3) |
---|
| 214 | def _weights(self, center, sigma, lb, ub): |
---|
| 215 | x = self._linspace(center, sigma, lb, ub) |
---|
| 216 | px = np.exp(-np.abs(x-center) / np.abs(sigma)) |
---|
| 217 | return x, px |
---|
[ce27e21] | 218 | |
---|
[ff7119b] | 219 | # dispersion name -> disperser lookup table. |
---|
[6cefbc9] | 220 | # Maintain order since this is used by sasview GUI to order the options in |
---|
| 221 | # the dispersion type combobox. |
---|
| 222 | MODELS = OrderedDict((d.type, d) for d in ( |
---|
| 223 | RectangleDispersion, |
---|
[75e4319] | 224 | UniformDispersion, |
---|
[6cefbc9] | 225 | ArrayDispersion, |
---|
| 226 | LogNormalDispersion, |
---|
| 227 | GaussianDispersion, |
---|
| 228 | SchulzDispersion, |
---|
[a5a12ca] | 229 | BoltzmannDispersion |
---|
[ce27e21] | 230 | )) |
---|
| 231 | |
---|
[f41027b] | 232 | SAS_WEIGHTS_PATH = "~/.sasview/weights" |
---|
[aa25fc7] | 233 | def load_weights(pattern=None): |
---|
| 234 | # type: (str) -> None |
---|
| 235 | """ |
---|
[f41027b] | 236 | Load dispersion distributions matching the given glob pattern |
---|
[aa25fc7] | 237 | """ |
---|
| 238 | import logging |
---|
| 239 | import os |
---|
| 240 | import os.path |
---|
| 241 | import glob |
---|
| 242 | import traceback |
---|
| 243 | from .custom import load_custom_kernel_module |
---|
| 244 | if pattern is None: |
---|
[f41027b] | 245 | path = os.environ.get("SAS_WEIGHTS_PATH", SAS_WEIGHTS_PATH) |
---|
| 246 | pattern = os.path.join(path, "*.py") |
---|
[aa25fc7] | 247 | for filename in sorted(glob.glob(os.path.expanduser(pattern))): |
---|
| 248 | try: |
---|
[910c0f4] | 249 | #print("loading weights from", filename) |
---|
[aa25fc7] | 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)) |
---|
[1780d59] | 254 | |
---|
| 255 | def get_weights(disperser, n, width, nsigmas, value, limits, relative): |
---|
[ff7119b] | 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 | |
---|
[6cefbc9] | 269 | *limits* is [lb, ub], the lower and upper bound on the possible values. |
---|
[ff7119b] | 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 | |
---|
[823e620] | 274 | Returns *(value, weight)*, where *value* and *weight* are vectors. |
---|
[ff7119b] | 275 | """ |
---|
[6cefbc9] | 276 | if disperser == "array": |
---|
[2d81cfe] | 277 | raise NotImplementedError("Don't handle arrays through get_weights;" |
---|
| 278 | " use values and weights directly") |
---|
[fa5fd8d] | 279 | cls = MODELS[disperser] |
---|
[1780d59] | 280 | obj = cls(n, width, nsigmas) |
---|
[823e620] | 281 | v, w = obj.get_weights(value, limits[0], limits[1], relative) |
---|
[a5a12ca] | 282 | return v, w/np.sum(w) |
---|
[6cefbc9] | 283 | |
---|
| 284 | |
---|
[a5a12ca] | 285 | def plot_weights(model_info, mesh): |
---|
| 286 | # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None |
---|
[6cefbc9] | 287 | """ |
---|
| 288 | Plot the weights returned by :func:`get_weights`. |
---|
| 289 | |
---|
[a5a12ca] | 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. |
---|
[6cefbc9] | 295 | """ |
---|
| 296 | import pylab |
---|
| 297 | |
---|
[2d81cfe] | 298 | if any(len(dispersity) > 1 for value, dispersity, weights in mesh): |
---|
[6cefbc9] | 299 | labels = [p.name for p in model_info.parameters.call_parameters] |
---|
[a5a12ca] | 300 | #pylab.interactive(True) |
---|
[6cefbc9] | 301 | pylab.figure() |
---|
[2d81cfe] | 302 | for (v, x, w), s in zip(mesh, labels): |
---|
[a5a12ca] | 303 | if len(x) > 1: |
---|
| 304 | pylab.plot(x, w, '-o', label=s) |
---|
[6cefbc9] | 305 | pylab.grid(True) |
---|
| 306 | pylab.legend() |
---|
[f1a8811] | 307 | #pylab.show() |
---|