Changeset ce27e21 in sasmodels for sasmodels/weights.py

Timestamp:

Aug 24, 2014 7:18:14 PM (10 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

1780d59

Parents:

14de349

Message:

first pass for sasview wrapper around opencl models

File:

• sasmodels/weights.py (modified) (1 diff)

sasmodels/weights.py

@@ +1 @@
+# TODO: include dispersion docs with the disperser models
+from math import sqrt
 import numpy as np
+from scipy.special import gammaln
 
-class GaussianDispersion(object):
-    def __init__(self, npts=35, width=0, nsigmas=3): #number want, percent deviation, #standard deviations from mean
-        self.type = 'gaussian'
-        self.npts = npts
-        self.width = width
-        self.nsigmas = nsigmas
+class Dispersion(object):
+    """
+    Base dispersion object.
+
+    Subclasses should define *_weights(center, sigma, lb, ub)*
+    which returns the x points and their corresponding weights.
+    """
+    type = "base disperser"
+    default = dict(npts=35, width=0, nsigmas=3)
+    def __init__(self, npts=None, width=None, nsigmas=None):
+        self.npts = self.default['npts'] if npts is None else npts
+        self.width = self.default['width'] if width is None else width
+        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas
 
     def get_pars(self):
-        return self.__dict__
+        pars = {'type': self.type}
+        pars.update(self.__dict__)
+        return pars
 
-    def get_weights(self, center, min, max, relative):
-        """ *center* is the center of the distribution
-        *min*,*max* are the min, max allowed values
-        *relative* is True if the width is relative to the center instead of absolute
-        For polydispersity use relative.  For orientation parameters use absolute."""
+    def set_weights(self, values, weights):
+        raise RuntimeError("set_weights is only available for ArrayDispersion")
+
+    def get_weights(self, center, lb, ub, relative):
+        """
+        Return the weights for the distribution.
+
+        *center* is the center of the distribution
+
+        *lb*,*ub* are the min and max allowed values
+
+        *relative* is True if the distribution width is proportional to the
+        center value instead of absolute.  For polydispersity use relative.
+        For orientation parameters use absolute.
+        """
         npts, width, nsigmas = self.npts, self.width, self.nsigmas
-        sigma = width * center if relative else width
+        sigma = self.width * center if relative else self.width
         if sigma == 0:
-            return np.array([center],'d'), np.array([1.], 'd')
-        x = center + np.linspace(-nsigmas * sigma, +nsigmas * sigma, npts)
-        x = x[(x >= min) & (x <= max)]
+            return np.array([center], 'd'), np.array([1.], 'd')
+        return self._weights(center, sigma, lb, ub)
+
+    def _weights(self, center, sigma, lb, ub):
+        """actual work of computing the weights"""
+        raise NotImplementedError
+
+    def _linspace(self, center, sigma, lb, ub):
+        """helper function to provide linear spaced weight points within range"""
+        npts, nsigmas = self.npts, self.nsigmas
+        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts)
+        x = x[(x >= lb) & (x <= ub)]
+        return x
+
+
+class GaussianDispersion(Dispersion):
+    type = "gaussian"
+    default = dict(npts=35, width=0, nsigmas=3)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, lb, ub)
         px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
         return x, px
 
+
+class RectangleDispersion(Dispersion):
+    type = "rectangle"
+    default = dict(npts=35, width=0, nsigmas=1.70325)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma*sqrt(3.0), lb, ub)
+        px = np.ones_like(x)
+        return x, px
+
+
+class LogNormalDispersion(Dispersion):
+    type = "lognormal"
+    default = dict(npts=80, width=0, nsigmas=8)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
+        px = np.exp(-0.5*(np.log(x)-center)**2)/sigma**2/(x*sigma)
+        return x, px
+
+
+class SchulzDispersion(Dispersion):
+    type = "schulz"
+    default = dict(npts=80, width=0, nsigmas=8)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8))
+        R= x/center
+        z = (center/sigma)**2
+        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
+        px = np.exp(arg)
+        return x, px
+
+
+class ArrayDispersion(Dispersion):
+    type = "array"
+    default = dict(npts=35, width=0, nsigmas=1)
+    def __init__(self, npts=None, width=None, nsigmas=None):
+        Dispersion.__init__(self, npts, width, nsigmas)
+        self.values = np.array([0.], 'd')
+        self.weights = np.array([1.], 'd')
+
+    def set_weights(self, values, weights):
+        self.values = np.ascontiguousarray(values, 'd')
+        self.weights = np.ascontiguousarray(weights, 'd')
+        self.npts = len(values)
+
+    def _weights(self, center, sigma, lb, ub):
+        # TODO: interpolate into the array dispersion using npts
+        x = center + self.values*sigma
+        idx = (x>=lb)&(x<=ub)
+        x = x[idx]
+        px = self.weights[idx]
+        return x, px
+
+
+models = dict((d.type,d) for d in (
+    GaussianDispersion, RectangleDispersion,
+    ArrayDispersion, SchulzDispersion, LogNormalDispersion
+))
+
+def get_weights(disperser, n, width, nsigma, value, limits, relative):
+    cls = models[disperser]
+    obj = cls(n, width, nsigma)
+    v,w = obj.get_weights(value, limits[0], limits[1], relative)
+    return v,w